Magnetorotational collapse of massive stellar cores to neutron stars: 
Simulations in full general relativity 
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We study magnetohydrodynamic (MHD) effects arising in the collapse of magnetized, rotating, 
massive stellar cores to proto-neutron stars (PNSs) . We perform axisymmetric numerical simulations 
in full general relativity with a hybrid equation of state. The formation and early evolution of a 
PNS are followed with a grid of 2500 x 2500 zones, which provides better resolution than in previous 
(Newtonian) studies. We confirm that significant differential rotation results even when the rotation 
of the progenitor is initially uniform. Consequently, the magnetic field is amplified both by magnetic 
winding and the magnetorotational instability (MRI). Even if the magnetic energy -Bem is much 
smaller than the rotational kinetic energy T rot at the time of PNS formation, the ratio E-EM/T TOt 
increases to 0.1-0.2 by the magnetic winding. Following PNS formation, MHD outflows lead to losses 
of rest mass, energy, and angular momentum from the system. The earliest outflow is produced 
primarily by the increasing magnetic stress caused by magnetic winding. The MRI amplifies the 
poloidal field and increases the magnetic stress, causing further angular momentum transport and 
helping to drive the outflow. After the magnetic field saturates, a nearly stationary, collimated 
magnetic field forms near the rotation axis and a Blandford-Payne type outflow develops along the 
field lines. These outflows remove angular momentum from the PNS at a rate given by J ~ jj-EemCb, 
where 77 is a constant of order ~ 0.1 and Cb is a typical ratio of poloidal to toroidal field strength. 
As a result, the rotation period quickly increases for a strongly magnetized PNS until the degree 
of differential rotation decreases. Our simulations suggest that rapidly rotating, magnetized PNSs 
may not give rise to rapidly rotating neutron stars. 

PACS numbers: 04.25.Dm, 04.30.-w, 04.40.Dg 



I. INTRODUCTION 

The explosion mechanism in core-collapse supernovae 
has been pursued for several decades, but the problem 
has not yet been solved. Neutrino-driven convection has 
been suggested as a way of rejuvenating the stalled shock, 
while rotation and magnetic fields have also been pro- 
posed as key driving mechanisms (see e.g., [1] for a re- 
view). Recently, acoustic power generated by accretion- 
triggered g-mode oscillations of the proto-neutron star 
(PNS) has been offered as an explanation for the explo- 
sion [2]. 

Even if rotation and magnetic fields are not the main 
mechanisms, they may still have an important influence 
on a supernova explosion, especially when a long-soft 
gamma-ray burst is produced. In addition, several ob- 
servations can be explained naturally by magnetic fields. 
For example, the rapid spindown of anomalous X-ray 
pulsars may result from the amplification of the star's 
magnetic field during core collapse [3, 4]. Soft gamma- 
ray repeaters [3] are likely to be neutron stars with very 
strong magnetic fields (so-called "magnetars" ) [5] . 

It has been speculated that the magnetic field strength 
could be amplified to ~ 10 16 G during stellar core col- 
lapse, and the resulting strong magnetic field could play 
a crucial role in the supernova explosion [6-9]. The col- 
lapse is generically nonhomologous, with the inner core 
collapsing faster than the outer core. Thus, even if the 
progenitor has a rigid rotation profile at the onset of the 



collapse, differential rotation naturally develops with an- 
gular velocity decreasing outward. In the presence of 
such differential rotation, the magnetic field is amplified 
both by magnetic winding [10-13] and the magnetoro- 
tational (magneto-shear) instability (MRI) [14-16]. The 
field strength may thus grow by many orders of magni- 
tude, even if it is initially small. The field amplification is 
likely to continue until the kinetic energy associated with 
the differential rotation Tdiff is converted to magnetic en- 
ergy E~em [13]. Typically, T^is is an appreciable fraction 
of the total rotational kinetic energy T rot . The amplified 
magnetic field results in a strong magnetic stress, which 
could blow off the matter in the vicinity of the PNS, con- 
verting magnetic energy back into matter kinetic energy 
and driving an outflow. The typical rotational kinetic 
energy of a PNS is approximately 

T rot = ^ Kl MR 2 n 2 

\0.3j\lAM Q ) 

x (ifjL) (S) ergs ' (1) 

where M, R, and fi(= 27r/P rot ) are the mass, radius, 
and angular velocity of the PNS, and kj ~0.3-0.6 de- 
notes the ratio of the moment of inertia to MR 2 and 
depends on the structure of the neutron star [17, 18]. 
Thus, if the rotation period of the PNS is shorter than 
~ 4 ms, and if a large fraction of the rotational kinetic 
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energy is converted to the kinetic energy of the outward 
matter flow via magnctohydrodynamic (MHD) processes, 
the supernova explosion may be significantly powered by 
magnetorotational effects. 

Numerical simulations of magnetorotational core col- 
lapse were pioneered by LeBlanc and Wilson [19] and by 
Symbalisty [20] in the 1970s and 1980s, respectively. In 
the past few years (e.g., [21-28]), this has become an ac- 
tive research topic in computational astrophysics. All of 
the simulations to date have been performed by assum- 
ing Newtonian gravitation or by including the general 
relativistic effects approximately [26]. Yamada and his 
collaborators have performed a variety of simulations in 
axial symmetry with simplified microphysics for a variety 
of rotational profiles and magnetic field strengths [21, 22]. 
They found that the initially poloidal magnetic field is 
primarily amplified by compression during infall and by 
magnetic winding following core collapse. They also 
found that the magnetic pressure, which is amplified dur- 
ing collapse, could drive a strong outflow along the ro- 
tation axis. Their simulations were done in the nonuni- 
form spherical polar coordinates (r, 9) with a grid size 
of at most (300, 50). At this resolution, they were not 
able to resolve the MM. Takiwaki et al. [23] performed 
simulations similar to those of Yamada et al., but with 
a realistic equation of state (EOS), and obtained essen- 
tially the same results. Kotake et al. [24] performed a 
simulation adopting a strong toroidal magnetic field and 
an extremely high degree of differential rotation as the 
initial condition and found that the toroidal magnetic 
field can power the explosion for this special initial con- 
dition. Obcrgaulingcr, Aloy, and Miillcr [25] performed 
similar simulations to those by Yamada et al. but with 
a slightly better grid resolution, (380, 60), in spherical 
polar coordinates. They indicated that not only mag- 
netic winding but also the MRI could play an important 
role in the supernova explosion and evolution of the PNS. 
However, the grid resolution they adopted was also not 
sufficient for resolving the regions in which the MRI oc- 
curs most strongly [29]. Ardeljan, Bisnovatyi-Kogan, and 
Moiseenko [27, 28] performed axisymmetric simulations 
with a Lagrangian scheme. In their work, a purely hydro- 
dynamic simulation was performed up to the formation 
of a PNS. Then they added a poloidal magnetic field to 
the PNS to study the MHD effects. They found that such 
effects, and, in particular the enhancement of the mag- 
netic field strength by magnetic winding, could result in 
buoyancy and trigger a supernova explosion. They also 
reported that an MRI-like instability played an impor- 
tant role. However, in their simulation, the instability is 
observed only after <~ 100P ro t, long after the toroidal field 
first becomes significant. But the fastest-growing mode 
for the (shear-type) MRI is amplified exponentially in a 
few rotation periods, irrespective of the strength of the 
toroidal fields [16]. Thus, the instability that they found 
is unlikely to be a shear- type instability [11]; rather it is 
likely to be the magneto-convective instability (see dis- 
cussion in Sec. VII). 



All of these previous simulations have provided a quali- 
tative picture of magnetorotational core collapse and su- 
pernova explosions. However, not all of the important 
magnetorotational effects have been incorporated in the 
simulations, mainly because of insufficient grid resolu- 
tion. Hence, some fundamental questions concerning the 
magnetorotational explosion scenario have not been an- 
swered. For example, it has been established that the 
magnetic field is amplified during the collapse and dur- 
ing the subsequent evolution of the PNS, but what de- 
termines the saturated field strength and how large is 
it? How important are the effects of the MRI? Rota- 
tional kinetic energy of the PNS can be converted to out- 
flow kinetic energy via MHD processes such as magnetic 
winding and the MRI. How efficiently is the rotational 
kinetic energy of the PNS converted to matter outflow 
energy? What is the rate of decrease of the rotational 
kinetic energy of the PNS and the corresponding rate of 
its period increase? After the magnetic field saturates, 
the magnetic configuration in the PNS will settle down 
to a stationary state. What is this final state? What 
mechanisms drive outflow? 

Previous work has been performed mainly in the frame- 
work of Newtonian gravitation. However, general rela- 
tivistic effects are not negligible and may have signifi- 
cant influence in the evolution of the PNS (particularly 
for a massive progenitor). Furthermore, after the col- 
lapse, shocks and outflows are often produced at rela- 
tivistic speed. From a technical point of view, general 
relativistic simulations have an additional advantage. In 
the Newtonian case, the Alfven velocity may exceed the 
speed of light (especially in the low-density region) , and 
hence, the Courant condition for the time step can be 
quite severe. On the other hand, the Alfven velocity is 
guaranteed to be smaller than the speed of light in gen- 
eral relativity, and so we do not have to deal with this 
unphysical complication. 

In this paper, we summarize results from axisymmet- 
ric simulations in full general relativity using two general 
relativistic magnetohydrodynamical (GRMHD) codes in- 
dependently developed recently [30-33]. One purpose is 
to demonstrate that multi-dimensional MHD simulations 
for core collapse are now feasible in full general rela- 
tivity, without approximation. The main purpose is to 
address the questions raised above by performing simu- 
lations with much higher grid resolution than those re- 
ported previously. We adopt a uniform grid in cylindrical 
coordinates so that we can achieve uniformly high reso- 
lution everywhere in our computation domain, and the 
MRI can be resolved wherever it occurs. We have also 
performed simulations using fisheye coordinates (see [34] 
and Appendix A), which provides high resolution in the 
central region while using relatively few of grid points. 
As in previous works [21, 22, 25], we employ a simpli- 
fied hybrid EOS since we want to focus on studying the 
MHD effects on the collapse and on the evolution of the 
magnetized PNS. 

The latest study on supernova progenitors [35] sug- 
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gests that magnetized, massive stars of solar metallicity 
may not be rotating as rapidly as the models we consider 
in this paper. It suggests that the typical value of P ro t 
at formation of a PNS is ~ 15 ms. The reason is that 
the magnetic field grows due to a dynamo mechanism 
during stellar evolution. As a result, angular momen- 
tum is transported outwards by magnetic braking, which 
decreases the rotation period in the central region of a 
massive star. If this scenario holds for all progenitors, 
the magnetorotational explosion scenario would not be 
effective since a substantial amount of rotational kinetic 
energy is required (P ro t 5 ms), as shown in Eq. (1). 
However, the rotational kinetic energy of the supernova 
progenitor depends strongly on model parameters of the 
dynamo theory and there is still a possibility of forming 
a rapidly rotating progenitor if the progenitor is mas- 
sive [35]. The magnetorotational scenario thus warrants 
detailed investigation. In this paper, we consider initial 
pre-collapse stars with substantial rotational kinetic en- 
ergies. 

The remainder of the paper is organized as follows. 
In Sec. II, we briefly review our mathematical formalism 
and numerical methods for our GRMHD simulations. In 
Sec. Ill, we provide a qualitative summary of the mag- 
netorotational effects that could play important roles in 
the evolution of a differentially rotating PNS. Sections IV 
and V describe our initial models and computational 
setup, respectively. In Sec. VI, we present our numerical 
results, focusing on the evolution of the magnetic fields 
and of the newly formed PNS. Finally, we summarize 
our conclusions in Section VII. Throughout this paper, 
we adopt geometrical units in which G = 1 = c, where 
G and c denote the gravitational constant and speed of 
light, respectively. Cartesian coordinates are denoted by 
x k = (x,y,z). The coordinates are oriented so that the 
rotation axis is along the z-dircction. We define the co- 
ordinate radius r = \/ x 2 + y 2 + z 2 , cylindrical radius 
vj = \J x 2 + y 2 , and azimuthal angle ip = tan~ 1 (y/x). 
Coordinate time is denoted by t. Greek indices /i, v, ■ ■ ■ 
denote spacetime components [x, y, z, and t), small Latin 
indices i,j, ■ ■ ■ denote spatial components (x,y, and z), 
and capital Latin indices /, J, ■ ■ ■ denote the poloidal 
components iw and z). 

II. FORMULATION 
A. Brief summary of methods 

The formulation and numerical scheme for the present 
GRMHD simulations are the same as those we reported 
in [30, 31], to which the reader may refer for details. 

The fundamental variables for the metric evolution 
are the three- metric 7^ and extrinsic curvature Kij. 
We adopt the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formalism [36-40] to evolve the metric. In 
this formalism, the evolution variables are the confor- 
mal factor ip = = 7 1 / 12 , the conformal 3-metric 



jij = e three auxiliary functions i* 1 , = S^ k djjik 

(or T l = — 7 u ,j)i the trace of the extrinsic curvature 
K, and the tracefree part of the extrinsic curvature 
Aij = e-^(K %J - jijK/3). Here, 7 = det( 7jJ ). The 
full metric <? M „ is related to the three-metric 7 M „ by 
Ifjtv = 9[j,v + n fl n u ^ where the future-directed, timclike 
unit vector n M normal to the time slice can be written in 
terms of the lapse a and shift (3 l as n M = a _1 (l, — 

The Einstein equations are solved in Cartesian coordi- 
nates. We employ the Cartoon method [41, 42] to impose 
axisymmetry. In addition, we also assume a reflection 
symmetry with respect to the equatorial plane and only 
evolve the region with x > and z > 0. We perform sim- 
ulations using a fixed uniform grid with size JVx3xiY 
in x — y — z which, covering a computational domain 
< x < L, < z < L, and -A < y < A. Here, N 
and L are constants and A = L/N. The variables in the 
y = ±A planes are computed from the quantities in the 
y = plane by imposing axisymmetry. 

The lapse a and shift (3 Z are gauge functions that have 
to be specified in order to evolve the metric. In the code 
of [31], an approximate maximal slice (AMS) condition 
K w is adopted following previous papers [37, 40, 43]. 
In this condition, a is determined by solving approxi- 
mately an elliptic-type equation. The shift is determined 
by the hyperbolic gauge condition as in [44, 45]. In the 
code of [30] , the lapse and shift are determined by hyper- 
bolic driver conditions as in [46]. 

The fundamental variables in ideal MHD are the rest- 
mass density p, specific internal energy e, pressure P, 
four- velocity u M , and magnetic field measured by an 
observer comoving with the fluid. The ideal MHD con- 
dition is written as u^F^ — 0, where F^ v is the elec- 
tromagnetic tensor. The tensor F^ and its dual in the 
ideal MHD approximation are given by 

F" v = t^ af3 u a b p , (2) 

F* v = ^e^ a pF af3 = - b v u^ (3) 

where e M „ Q /3 is the Levi-Civita tensor. The magnetic field 
measured by a normal observer n M is given by 

= n v F*'^ = aiu'b^ - b'u"), (4) 

From the definition of B^ and the antisymmetry of F^ u 
we have n^B^ = = 5*. The relations between 6 P and 
B l are 



a 

b l = -L(b* + Biu 1 u l ). (6) 

Using these variables, the energy-momentum tensor is 
written as 

rp rTiFlllid I rpEM (rj\ 

± iiv — J- ^ u i-'/iv 1 \' ) 

where TjJ uld and r^ M denote the fluid and electromag- 
netic pieces of the stress-energy tensor. They are given 
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by 



T™ d = phu fl u„ + Pg fiV , 

T EM p p a _ 1 F aF alS 



(8) 



(9) 



Here, ft. = 1 + e + P/p is the specific enthalpy and b 2 = 
b^b^. Hence, the total stress-energy tensor becomes 

= (ph + b 2 )u^u v + (p + h -^j 9llv - b^b u . (10) 

The quantity b 2 is often referred to as the magnetic en- 
ergy density, and b 2 /2 = P mag as the magnetic pressure. 

In our numerical implementation of the GRMHD and 
magnetic induction equations, we evolve the weighted 
density , weighted momentum density Si , weighted en- 
ergy density So, and weighted magnetic field B % . They 
are defined as 



P* = — y/ypn^, 

So = y/iT^rfn", 
B l = 



(11) 
(12) 

(13) 

(14) 



During the evolution, we also need the three- velocity v % = 
u 1 ju l . 

The GRMHD and induction equations are written in 
conservative form for variables p», Si, So, and B % and 
evolved using a high-resolution shock-capturing scheme 
(HRSC). In the code of [31], the high-resolution central 
(HRC) [47, 48] scheme is employed. In this approach, 
the transport terms such as • •) are computed by 
a Kurganov-Tadmor scheme with a third-order (piece- 
wise parabolic) spatial interpolation. It has been demon- 
strated [48] that the results obtained using this scheme 
approximately agree with those using an HRSC scheme 
based on Roe- type reconstruction for the fluxes [42, 52]. 
In the code of [30], the transport terms are evaluated by 
the HLL (Hartcn, Lax and van-Leer) flux formula [49]. 
The cell interface data reconstruction is done by the 
monotonized central (MC) scheme [50]. It has also been 
demonstrated that the performance of the HLL scheme is 
as good as the Roe- type and HRC schemes in many MHD 
simulations [51]. The magnetic field B l has to satisfy the 
no monopole constraint diB % = 0. In the code of [31], 
this magnetic constraint is imposed by using the con- 
strained transport scheme first developed by Evans and 
Hawley [31, 53]. In the code of [30], the flux-interpolated 
constrained transport (flux-CT) scheme [54] is employed. 
Both constrained transport schemes guarantee that no 
magnetic monopoles will be created in the computation 
grid during the numerical evolution. At each timestep, 
the primitive variables (p,P,v l ) must be computed from 
the evolution variables (p*,So,Si). This is done by nu- 
merically solving the algebraic equations (11)— (13) to- 
gether with the EOS P = P(p,e). As in many hydro- 
dynamic simulations in astrophysics, we add a tenuous 



"atmosphere" that covers the computational grid out- 
side the star. The atmospheric rest-mass density is set 
to w 10~ 10 p c (0) (w 1 g/cm 3 ), where p c (Q) is the initial 
central density of the star. 

The codes used here have been tested in relativis- 
tic MHD simulations, including MHD shocks, nonlin- 
ear MHD wave propagation, magnetized Bondi accretion, 
MHD waves induced by linear gravitational waves, and 
magnetized accretion onto a neutron star. Furthermore, 
we have used both codes to perform simulations of the 
evolution of magnetized, differentially rotating, relativis- 
tic, hypermassive neutron stars, and obtain good agree- 
ment [32, 33]. 



B. Equations of state 

A parametric, hybrid EOS is adopted following Miiller 
and his collaborators [25, 55, 56]. In this EOS, the pres- 
sure consists of the sum of a cold part and a thermal 
part: 



P(p,e) = P P (p) + P th (p,e). 



(15) 



The cold part of the pressure Pp depends only on the 
density. In this paper, we choose the following form of 



■ K x p Tl , p < p nuc , 



(16) 



where K\, K 2 , V\ and T2 are constants, and p nuc w 
2 x 10 14 g/cm 3 is nuclear density. We set K 2 = Kip^~ F2 
to make Pp continuous at p = p nuc . The adiabatic in- 
dices are chosen as Ti = 1.3 or 1.32 and F 2 =2.5 or 2.75. 
Following [55, 56], the value of K\ is chosen to be 5 x 10 14 
in cgs units. With this choice, the cold part of the EOS 
for p < p nuc is approximately given by relativistic de- 
generate electron pressure. This simplified cold EOS is 
designed to mimic a more complicated cold nuclear EOS. 
Using the first law of thermodynamics at zero tempera- 
ture, we obtain the specific internal energy ep associated 
with the cold part of the pressure Pp: 



e P (p) = - J PAp)d(^j 



= < 



Ti-1' 

l r 2 - 1' 



P < Pr, 



r ., , (r 2 - ro^pii 



r,-i 



(17) 



+ 



(Tx - i)(r 2 - l) 



P > Pn 



The thermal part of the pressure P t h plays an important 
role when shocks occur. We adopt a simple form for P t h- 



P 



th 



(r t h - l)/0£th, 



(18) 



where £th = e — £p is the thermal specific internal energy. 
The value of T t h determines the efficiency of converting 
kinetic energy to thermal energy at shocks. We set r t h = 
Ti to conservatively account for shock heating. 
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C. Diagnostics 

We monitor the total baryon rest mass M», ADM 
(Arnowitt-Deser-Misner) mass M, and angular momen- 
tum J, which are computed as in [33, 42]. We also com- 
pute the internal energy E- lnt , thermal internal energy 
T-hcat , rotational kinetic energy T rot , gravitational poten- 
tial energy W, and electromagnetic energy -Eem using the 
formulae given in [33]. 

In a stationary, axisymmetric spacctimc, the energy E 
and angular momentum J are conserved, where 

E = j a^ffT^x, (19) 
J = J a^jT^cPx. (20) 

For a nearly stationary spacetime, which is achieved af- 
ter the formation of the PNS, E and J are approximately 
conserved. We can then define the fluxes of rest mass, 
energy, and angular momentum across a sphere of coor- 
dinate radius r by 

F M (r) = <f dAp*v r , (21) 

J r— const 

F E (r) =-<f dAa^T r t , (22) 

J r— const 

Fj(r) = <f dA a ^T r v , (23) 

J r— const 

where dA = r 2 sin 6d6d(j>. The energy and angular mo- 
mentum fluxes associated with the electromagnetic fields 
are defined as 

F E .VM{r) =-<f dAa^T EMr t , (24) 

J r— const 

fj,EM(r) - / dAa^T™ r v . (25) 

J r— const 

The total energy flux Fe is very close to the rest-mass 
flux F M since F E is primarily composed of the rest-mass 
energy flow. Thus, we define another energy flux by sub- 
tracting the rest-mass flow: F e = F E — Fm- Wc note 
that F e contains kinetic, thermal, electromagnetic, and 
gravitational potential energy fluxes. If F e > at suf- 
ficiently large radius, an unbound outflow (overcoming 
gravitational binding energy) is present. 

It should be noted that even in a stationary spacetime, 
one can choose gauges such that d/dt is not a Killing 
vector. In this case, the energy E defined above is not 
conserved. Thus, these are physical meaningful fluxes 
only in certain gauges. In our evolution, we find that 
after the formation of the PNS, the system relaxes to a 
stationary state and the metric does not change signifi- 
cantly with coordinate time t, suggesting that d/dt is an 
approximate Killing vector in our simulations. Hence we 
use the above formulae to compute the fluxes. 



D. Gravitational waveforms in terms of quadrupole 
formula 

We compute gravitational waves in terms of the 
quadrupole formula given by [58, 59]: 

"«=^ p /(^> < 26 » 

where 7^ and P i 3 = Sij — hihj (hi = x l /r) denote 
the tracefree quadrupole moment and the TT projection 
tensor, respectively. From this expression, the +- mode 
of quadrupole gravitational waves in an axisymmetric 
spacetime can be written as 

hf ad = J **fc°tW«(*rct) gin 2 ^ (2?) 

r 

where 7y denotes the quadrupole moment, is its sec- 
ond time derivative, 9 denotes the angle between the ro- 
tation axis and the direction of observation, and t let is 
retarded time i rot ~ t — r. In this paper, we charac- 
terize the gravitational waves by the variable A2(t) = 
I X x{tret) — Izz(trct), which has dimensions of length. In 
terms of A 2 , the observed gravitational- wave strain is 
given by 

In spacetimcs with strong gravitational fields, there is 
no unique definition for the quadrupole moment. Here, 
we choose for simplicity 

hj = j p*x i x j d 3 x. (29) 

Using the continuity equation, we compute the first time 
derivative as 

4 ■= J P*{v l x J +x l v j )d 3 x. (30) 

We then compute 4 by finite differencing the numerical 
values of 7^ in time. 

Because of the ambiguity in the definition of the 
quadrupole moment in a spacetime with strong gravita- 
tional fields, one may choose alternative expressions for 
Iij . Hence, the gravitational waveforms determined from 
the quadrupole formula depend on the chosen definition 
of Iij. In addition, they depend on the gauge conditions, 
since the coordinates r and t appeared in 7^ and 7jj are 
gauge dependent. We have calibrated our quadrupole for- 
mula in [58], and found that the magnitude of the error 
in the amplitude of gravitational waves is of order M / R, 
where M and R denote the typical mass and radius of the 
emitter. On the other hand, the phase of the waveforms 
is quite accurate. Thus, the wave amplitudes shown in 
this paper are not accurate to better than ~10%, but 
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the computed radiation possesses the correct qualitative 
features. 

We also note that in our formula, the contribution from 
is neglected. This is justified in the present treat- 
ment, since the contribution from the electromagnetic 
part is about 10% of the matter part, and hence, the er- 
ror in neglecting the former is as large as the error of our 
quadrupole formula. 



III. MAGNETOROTATIONAL EFFECTS 

In this section, we summarize processes that play an 
important role in magnetorotational core collapse and the 
subsequent evolution of the PNS. 



A. Compression and magnetic winding during 
collapse 

The magnetic field is amplified by compression and 
magnetic winding during stellar collapse. This can be 
understood by considering the magnetic induction equa- 
tion in a perfectly conducting (ideal MHD) plasma: 

d t B l + d j {v j B i - v i B j ) = 0. (31) 

Combining Eq. (31) with the continuity equation 



Hence for an observer comoving with the fluid, 



dp* i P* o / I\ n 

— - H di{wv ) = 0, 

at zd 

we rewrite the induction equation as 
dH 1 



dt 
dW 
dt 



= 'H J djv I , 



where W = B l /p* 7 = O, and 

d_ = d_ j_d_ 

dt~ dt +V dx 1 ' 



(32) 

(33) 
(34) 

(35) 



For a pre-collapse star having moderate angular mo- 
mentum, the core collapses in an approximately spherical 
and homologous manner. Hence v r oc r, which implies 
w ro oc zd, v z oc z, and s=a \v z \. It follows from the 
continuity equation that 



dp, „ w ro 
— w -3p* — . 
dt zd 



(36) 



Note that t> ro and v z are negative. Combining the above 
equation with Eq. (33), we obtain 



dt zd 
dt zu 



rflnp* 
3 df~ 
H z d In p* 
~3 dt~' 



(37) 
(38) 



n 1 



oc p* 



1/3 



-1/3 



(39) 



which means that the poloidal field B 1 oc p 2 / 3 during 
core collapse. 

While the collapse proceeds roughly homologously in 
the bulk of the star, this is not the case in the outer lay- 
ers. Thus, significant differential rotation develops in the 
outer layers. We see from Eq. (34) that toroidal mag- 
netic fields are amplified when there is differential ro- 
tation along the poloidal field lines (magnetic winding). 
The toroidal field is thus expected to grow during the 
collapse. However, since the growth depends on the non- 
homologous nature of the collapse in the outer layers, a 
simple prediction for the dependence of Ti v on p [as in 
Eq. (39)] is not available. 

To obtain a rough idea of how the toroidal field evolves 
during the collapse, we proceed as follows. Assume that 
f2 in the differentially rotating outer regions is described 
by tu~ p , where p is constant with time. This will not 
apply in regions where the homology is seriously violated. 
However, the true behavior of in the region of interest is 
not known analytically, and we are only seeking a rough 
estimate. Equation (34) gives 



dH v 



-p- 



(40) 

dt zu K J 

The evolution of H v measured by the comoving observer 
is given by 



I 



dtp- 



zu 



(41) 



For a comoving observer, zd^ 1 oc p 1 / 3 and hence Ti m /zd 
is approximately constant according to Eq. (39). (This 
again assumes homologous collapse, which does not 
strictly hold in the regions of interest for the toroidal 
field growth.) When the magnetic field is weak and the 
spacetime is axisymmetric, the specific angular momen- 
tum of a fluid particle j — hu v w Qzu 2 is approximately 
conserved, which implies oc zd~ 2 oc p 2//3 . We then have 



H v oc 



Jndttx J p 3/2 dt. (42) 

Next, we obtain an approximate relation between dt and 
dp from the following Newtonian analysis. Since the col- 
lapse is approximately spherical and homologous, conser- 
vation of energy implies 



1 



v 2 _ m 0) 



m(ro 
ro 



(43) 



where ro is the radius of the fluid particle at the onset 
of collapse (t = 0). The mass interior to the radius r, 
m(r) — m(r ) — (4-7r/3)rQp is independent of time dur- 
ing homologous collapse. Here po is the density of the 
fluid particle at t — 0. Setting v = —dr/dt, we have 



dr 
~dl 



' m(ro) 



1 



1 

ro 



(44) 
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Combining the above equation with p = p (r /r) 3 , we 
obtain 

dtoc^V* ^ . . (45) 

(p/po) 4 ' 3 y/(p/Po) 1/3 1 

Substituting Eq. (45) into Eq. (42), we have H v oc 

V (p/po) 1 ^ 3 — 1j an d hence we provide the following es- 
timate of the toroidal field during the collapse: 

When p 3> po , the above equation simplifies to 

H v o:p 1/6 and B v cx p 7 ' & . (47) 

In practice, we find that B v cx p 9 , where q is slightly less 
than f. However, with the assumptions that went into 
Eq. (47), only rough agreement is to be expected. 



period, and wc have assumed a Keplerian angular veloc- 
ity profile f2 cx n7~ 3 / 2 . The growth of B T is expected to 
deviate from this linear relation when the magnetic ten- 
sion is large enough to change the angular velocity profile 
of the fluid (magnetic braking) . Magnetic braking trans- 
ports angular momentum on the Alfven timescale [10- 
13]: 




where R is the characteristic radius of the PNS and va ~ 
\B\/p x / 2 is the Alfven speed. So Eq. (51) holds for t p < 
t<t p + t A . 



C. Magnetorotational (shear) instability 



B. PNS magnetic winding 

When the central density of the collapsing star exceeds 
the nuclear density, the core bounces due to the stiff nu- 
clear EOS. The star quickly settles down to a quasista- 
tionary state and a PNS is formed. In this quasistation- 
ary state, the toroidal magnetic field grows linearly with 
time. This can be seen from the induction equation (31). 
If the magnetic field is weak and has a negligible back- 
reaction on the fluid, the velocities will remain constant 
with time. In cylindrical coordinates, we have (assuming 
axial symmetry) 



dtB 1 « 0, 

d t B v w —d^wQB 1 ), 



(48) 
(49) 



where we have assumed that <C zufl and \v z \ <C tail 
when the PNS is in quasiequilibrium. Then, 



d t B v w B'd^, 



(50) 



where we used the no-monopole constraint \di(wB I ) = 
in axisymmctry] . During the early phase of the PNS 
evolution, Eq. (50) indicates that the toroidal component 
of the field B T (— wB v ) grows linearly according to 

\B T (t;w,z)\ 

w \B T {t p ; vj,z)\ + tmlB^tp] zu, z)difl(t p ; zu, z)\ 
w \B T (t p ;w,z)\ + — \B™(t p ;m,z)\ 

"rot 



\B T (t p] w,z)\ 



+ 10 



15 



t 



100 ms 7 V 1 ms 



\B™{t p -zu lZ ) 
10 12 G 



qsi) 



where t p is the time at which the PNS first settles down to 
a quasiequilibrium state, P ro t denotes the local rotation 



The MRI is present in a weakly magnetized, rotat- 
ing fluid wherever <9 ra f2 < [14-16]. When the instabil- 
ity reaches the nonlinear regime, the distortions in the 
magnetic field lines and velocity field lead to turbulence. 
To estimate the growth time scale £mri and the wave- 
length of the fastest-growing mode A max , we can use a 
Newtonian local linear analysis given in [16]. Linearizing 
the MHD equations for a local patch of a rotating fluid 



on 



and imposing a plane-wave dependence (e 
the perturbations leads to the dispersion relation given 
in Eq. (125) of [16]. Specializing this equation for a 
constant-entropy medium leads to the dispersion relation 

w 4 - [2(k-v A ) 2 + K V 

+ (k • v A ) 2 [(k • V4) 2 + k 2 - 4fl 2 ] = , (53) 

where va = JS/yfp is the (Newtonian) Alfven velocity 
vector and k is the epicyclic frequency of Newtonian the- 
ory: 



1 d{vu 4 n 2 ) 



ZU' J 



dvu 



(54) 



Equation (53) is modified for a medium of inhomoge- 
neous entropy [16]. In this section, we neglect any en- 
tropy gradients and focus on the effects of shear. Wc 
further simplify the analysis by considering only the ver- 
tical modes (k = ke z ) since these are likely to be the 
dominant modes. The value of to 2 can be found by solv- 
ing Eq. (53) and then minimized to obtain the frequency 
of the fastest-growing mode, w max : 



2 _ 1 / dfl 
This maximum growth rate corresponds to 

-4 



(kv A ) 



z \2 
max 



16ft 2 



(55) 



(56) 
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The growth time (e-folding time) and wavelength of the 
fastest-growing mode are then given by 



*MRI = l/(iWmax) = 2 \dQ,/dlllw\ 



2tt 



2TTV* 



' 2il 



-1/2 



(57) 
(58) 



For a Keplerian angular velocity distribution Q cx w 3 / 2 , 
we have 

<MM = 4 " ° mS ( lOQols-0 ' (M) 
8ttw^ 



A max — 



i5fi 



2.1 km 







1000 rad s" 1 



10 14 G y VlO 13 g/cm 3 



-1/2 



.(60) 



If £F is comparable to the field strength of a canonical 
pulsar (B z ~ 10 12 G), A max is much smaller than the 
typical radius R of a PNS. Since A max cx va, larger mag- 
netic fields will result in longer MRI wavelengths. When 
Amax <5 100 km, the MRI will be suppressed since the 
unstable perturbations will no longer fit inside the re- 
gion with a high degree of differential rotation. Hence 
the MRI is regarded as a weak-field instability. In this 
paper, the initial magnetic field strength is chosen so that 
the field strength in the resulting PNS is ~ 10 14 G, giving 
Amax ~ a few km. 

Unlike A ma x, ^mri does not depend on the magnetic 
field strength but on the angular velocity profile. The 
Newtonian local analysis indicates that the MRI always 
grows on the timescale of a rotation period for a con- 
figuration with d In fl/d In w ~ — 1. Hence, the MRI is 
expected to play an important role for a PNS, especially 
near its surface, which is in general differentially rotating. 
The resulting strong magnetic fields and turbulence tend 
to transport angular momentum from the rapidly rotat- 
ing inner region of the PNS to the more slowly rotating 
outer layers. This causes the inner part to contract and 
the outer layers to expand. We note, however, that once 
the magnetic field strength is saturated, the resulting an- 
gular momentum transport and matter outflow will be 
governed by the turbulence and is thus expected to oc- 
cur on a time scale longer than £mri (e.g., a turbulent 
transport time scale). 



D. MHD outflow by the magneto-spring effect 

After the formation of the PNS, MHD outflows may 
transport angular momentum outward (causing a spin- 
down of the PNS) and may power a supernova explo- 
sion [6, 7]. In the early stage of the evolution of the PNS, 
the MHD outflow may be driven by the toroidal magnetic 



field, which grows linearly due to magnetic winding ac- 
cording to Eq. (51). A helical magnetic field forms as a 
result. The rate of change of the magnetic energy asso- 
ciated with the toroidal field is 



d(B T ) 2 
dt 

2B T B™zud zo n w -3B T B™n, 



E 



EM 



2B T B T 



(61) 



where we have assumed a Keplerian angular velocity pro- 
file cx ti7~ 3 / 2 . We note that B T B™ should be negative 
[see Eqs. (34) and (50)]. Because of the growing mag- 
netic pressure and hoop stress, material is lifted from the 
PNS along the rotation axis, producing a tower-like MHD 
outflow (see e.g., [63-65]). 

If a significant amount of energy gained from magnetic 
winding is used to drive the MHD outflow, the rate of 
energy loss from the PNS is approximately 



E ~ E] 



EM 



B z 



n. 



(62) 



where we have assumed that the magnetic energy is dom- 
inated by the toroidal field B T . This process is likely to 
be important only in the early stage before the toroidal 
field saturates. 



E. MHD outflow by magneto-centrifugal effect 

Another type of MHD outflow may occur in a differ- 
entially rotating object penetrated both by poloidal and 
toroidal field lines. Blandford and Payne showed that 
when a stationary disk is penetrated by open magnetic 
fields and when a mechanism for matter ejection operates 
at the disk surface, a stationary magneto-centrifugal out- 
flow is launched [7, 66] (see also relevant computational 
work, e.g., [67, 68]). 

We first consider the magneto-centrifugal effect oper- 
ating on an accretion disk. To estimate the order of 
magnitude of the energy flux in the outflow, we evalu- 
ate the Poynting flux approximately. We assume that 
the disk rotates with an angular velocity fid and that 
the poloidal velocity is smaller than the rotational ve- 
locity The Poynting flux along the poloidal field line 
is <~ \E x B\ ~ Rd£ldB p B T , where Rd is the charac- 
teristic radius of the disk, B p ~ _B ro is the magnitude 
of the poloidal magnetic field, and the electric field E 
is calculated by the (Newtonian) ideal MHD condition 
\E\ = \v x B\ ~ RdfldB p . The energy generation rate is 
then given by 



-Epoyn = (Poynting flux) x (area) 
~ Vd B p B T R d n d 



cx E EM C B (-^- 



(63) 



where r\d is a constant of order unity, Hd is the thickness 
of the disk, and Cb denotes a typical ratio of |-B P | to 
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\B T \. Thus, the rate of energy loss of the disk is similar 
to that of the magneto-spring mechanism. 

In the Blandford and Payne model [66], the Poynting 
flux is accompanied by a kinetic energy flux of compa- 
rable amplitude from matter outflow in the vicinity of 
the disk. Hence, the total rate of energy loss from the 
disk can be written in the same form of the last line of 
Eq. (63). 

The same qualitative analysis may be applied to a PNS 
by setting Hd ~ Rd, giving the energy loss rate as 



E = t]*EemCbQ, 



(64) 



where fi is the typical angular velocity of the PNS and 
77* is a constant of order unity. 

Assuming that the main source for the MHD outflow 
(cither by the magneto-spring or magneto-centrifugal 
mechanism) is the rotational kinetic energy of the star, E 
is related to the angular momentum flux J approximately 
by 



and hence, 



E w fij, 



J w t/^EemCb- 



(65) 



(66) 



For an MHD outflow, E and J should be approxi- 
mately equal to the fluxes integrated over a closed surface 
far away from the source. Thus, in this paper, we cal- 
culate Fm, F e , and Fj to obtain the rate of loss of the 
energy and angular momentum of the PNS. 

In the model of Blandford and Payne [66] , a priori mass 
ejection from the central object at a launching velocity 
as large as the rotational velocity is required to drive 
the MHD outflow. Thompson et al. [8] suggest that a 
neutrino wind could be the engine of the initial velocity. 
Here we anticipate a purely MHD source. In the present 
context, the magnetic field in the vicinity of the PNS 
is amplified by magnetic winding and the MRI. The re- 
sulting Poynting flux propagates outward and injects the 
energy into the Blandford-Payne wind. 



IV. PRE-COLLAPSE STELLAR MODELS 

The pre-collapse stars are modeled as rotating T = 
4/3 polytropes. Following [56, 59], the central density is 
chosen to be p c — 10 10 g/cm 3 . The EOS is given by 



(67) 



where K is set to be 5 x 10 14 in cgs units. This EOS cor- 
responds to the degenerate pressure of ultra-relativistic 
electrons [57]. We adopt a commonly used angular ve- 
locity profile [60, 61]: 



u u. 



(68) 



where O c denotes the angular velocity along the rotation 
axis, and zud is a constant. In the Newtonian limit, this 
rotation law reduces to the so-called 'j-constant' law: 



= a. 



' W -\- ZUj 



(69) 



Hence the parameter Wd controls the degree of differen- 
tial rotation of the star. In this paper, we choose rigidly 
rotating cases (zud — > 00) and a moderately differen- 
tially rotating case with A = Tu d /R = 0.5, where R is 
the equatorial radius. In the rigidly rotating cases, we 
choose the ratios of polar to equatorial radii, z s /R, to 
be 0.667 and 0.883. With z s /R = 0.667, the angular ve- 
locity at the equatorial surface is approximately equal to 
the Keplerian velocity (i.e., a uniformly rotating star at 
the mass-shedding limit). These two models with rapid 
and moderate rotation are referred to as models A and 
B, respectively. For the differentially rotating case, we 
choose a model with the ratio of the rotational kinetic 
energy T ro t to the gravitational potential energy W(> 0) 
of Rj 0.009, which is approximately the same as that of 
model A. In Table I we summarize the parameters of our 
models. 

To induce collapse, we fix the density distribution of 
the star but recalculate the pressure P and specific inter- 
nal energy s using the hybrid EOS [Eqs. (15)— (18)] . We 
set K~i = K p^ 3 ~ Fl , where po = 1 g/cm 3 . In a previous 
study of core collapse [59] , the pressure is computed by 
P = Kip Tl following [56]. The specific internal energy is 
then fixed by the hybrid EOS to be (note that we choose 

r th = r x ) 



c = 



(70) 



However, this choice reduces e by ^ 50% from the equi- 
librium value for Ti = 1.3, and the collapse is strongly 
accelerated in the early times. As a result, the collapse 
is strongly non-homologous and hence not very realistic. 
Therefore, in this paper, we do not decrease e, i.e., we 
set it according to 



e = 3K a p^ 3 . 



(71) 



The pressure is then determined by the hybrid EOS to 
be 



3^o (ri 



r 



n 4 /3 



(72) 



which is reduced uniformly from its initial value [cf. 
Eq. (67)] by a factor of 1 - 3(r x - 1) (i.e., 10% for 
Ti = 1.3). With this choice, the collapse is more uni- 
form, although the collapse time is longer. 

Next, we add a magnetic field to the stars. Since 
the initial profile of the magnetic field in the core col- 
lapse progenitor is not known, we follow our previous 
papers [32, 33] and add a dipole-like poloidal magnetic 
field to the pre-collapse model by introducing a vector 
potential of the following form: 



A v = AbTu 2 m&x{p 



l/n b 



pILTA- 



(73) 
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TABLE I: Central density p c , baryon rest mass M», ADM mass M, equatorial radius R, ratio of the rotational kinetic energy 
to the potential energy T Iot /W , non-dimensional angular momentum parameter J/M 2 , central value of the lapse function a c , 
angular velocity at the rotation axis Q c , ratio of the angular velocity at the rotation axis to that at the equatorial surface 
n c /Q s , ratio of polar to equatorial radii z s /R, and A — vj d /R of the pre-collapse stars. 



Model 


p c (g/cm 3 ) 


M*(M G ) 


M(M Q ) 


R (km) 


Trot/W 


J/M 2 


OLc 


fi c (rad/s) 


Qc/fls 


z s /R A 


A 


1.00 x 10 10 


1.503 


1.503 


2267 


8.9 x 10" 3 


1.235 


0.994 


4.11 


1.0 


0.667 oo 


B 


1.00 x 10 10 


1.486 


1.486 


1445 


4.9 x 10~ 3 


0.909 


0.994 


3.12 


1.0 


0.883 oo 


C 


1.00 x 10 10 


1.504 


1.504 


1628 


9.0 x 10" 3 


1.187 


0.994 


6.31 


5.0 


0.917 0.5 



Here the cutoff density p cut is chosen as 10~ 4 p c = 
10 6 g/cm 3 . (Simulations have been performed for a dif- 
ferent cutoff density, 10 8 g/cm 3 , and we find that the re- 
sults depend only weakly on this value.) The parameter 
Ab determines the initial strength of the magnetic field 
(see below), whereas rib in Eq. (73) determines the ini- 
tial profile of the magnetic field lines. We choose rib — 1 
and 4. In Fig. 1, we show the density contour curves and 
the magnetic field lines (which coincide with contours of 
A v in axisymmetry) for model A with rib = I (left) and 
with rib = 4 (right). This figure shows that for the larger 
value of rib, the maximum of the magnetic field strength 
B max is located at a larger cylindrical radius. This re- 
sults in a larger magnetic energy and a larger value of 
(•Pmag/-P)max for a given value of S max . 

We choose Ab such that the z-component of the mag- 
netic field is 7x 10 12 -4x 10 13 G. We summarize in Table II 
the maximum value of \B Z \, the ratio of the magnetic en- 
ergy to the rotational kinetic energy £?em /^rot ; arid the 
maximum value of Pm ag /P for the models we consider. 
With these values of Ab, the strength of the poloidal 
magnetic field of the resulting PNSs can be <; 10 15 G 
and the value of A max for the fastest-growing mode of 
the MRI becomes ^ 1 km, which can be resolved with 
our computational resources. Furthermore, the Alfven 
timescale in the PNSs is short enough (<J 30 ms) to see 
the whole magnetic braking process with a reasonable 
computational cost. 

Although the initial value of the magnetic field 
strength is large (much larger than that of a presupernova 
stellar model [35] in which the magnetic field is assumed 
to be amplified primarily by a dynamo process), the mag- 
netic energy is only ~ 0.1-5% of the rotational kinetic en- 
ergy, and is an even smaller fraction of the gravitational 
potential energy and the internal energy. Even if the sim- 
ulation starts with a smaller magnetic field strength, the 
fields are amplified globally by magnetic winding and lo- 
cally by the MRI until they saturate (see Sec. VI). Thus, 
the qualitative features of the evolution of the PNS may 
not depend strongly on the initial magnetic field strength. 

According to the observed number ratio of magnetars 
to canonical pulsars, the lower limit of the Galactic mag- 
netar birth rate exceeds 0.01 /year, which is as large as 
that of the canonical pulsars [3]. Namely, neutron stars 
with magnetic field strength ~ 10 15 G are not rare in the 
Universe. The origin of magnetars is not clear. A simple 
hypothesis is that the stellar progenitor of a magnetar 
has a strong magnetic field. The simulations in this pa- 



TABLE II: Parameters of the initial magnetic fields. 



Model 


n b 


\B Z |max(G) 


-EEIVl/Trot 


(Pm, S /P) 


max 


AO 

















Al 


1 


7.2 x 10 12 


1.9 x 10" 


-3 


2.1 x 10 


-4 


A2 


1 


1.8 x 10 13 


1.2 x 10" 


-2 


1.3 x 10 


-3 


A3 


1 


3.6 x 10 13 


4.8 x 10" 


-2 


5.2 x 10 


-3 


A4 


4 


7.2 x 10 12 


2.3 x 10" 


-2 


1.7 




B 


1 


1.8 x 10 13 


2.1 x 10" 


-2 


1.3 x 10 


-3 


C 


1 


7.2 x 10 12 


1.9 x 10" 


-3 


2.1 x 10" 


-4 



per provide such a model for the formation of magnetars. 



V. GRID SETUP 

During the collapse, the central density increases from 
10 10 g/cm 3 to - 6 x 10 14 g/cm 3 . This implies that the 
characteristic length scale of the system varies by a fac- 
tor of ~ 100. One of the computational challenges in a 
stellar core collapse simulation is maintaining numerical 
accuracy despite the significant change in the character- 
istic length scale. In the early phase of the collapse (infall 
phase; see Sec. IV A), which proceeds in a nearly homolo- 
gous manner, we may follow the collapse with a relatively 
small number of grid points, and successively move the 
outer boundary inward (while keeping the same number 
of grid points) to increase resolution. As the collapse 
proceeds, the central region shrinks more rapidly than 
the outer region, and hence a higher grid resolution is 
necessary to accurately follow such rapid collapse in the 
central region. On the other hand, the location of the 
outer boundaries cannot be changed by a large factor 
because this would discard matter in the outer envelope. 

To follow the collapse accurately and save CPU time, 
we adopt the regridding technique described in [59, 62]. 
Regridding is carried out whenever the characteristic ra- 
dius of the collapsing star decreases by a factor of 2-3. At 
each regridding, the grid spacing is decreased by a factor 
of 2. All the quantities in the new grid are calculated us- 
ing cubic interpolation. To avoid discarding the matter 
in the outer region, we also increase the grid number at 
regridding, ensuring that the discarded baryon rest mass 
is less than 3% of the total. 

Specifically, the values of N and L (see Sec. II A) 
in the present work are chosen in the following man- 
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FIG. 1: Initial density contour curves and magnetic field lines for model A with n& = 1 (left) and with rib = 4 (right). The 
density contour curves (thick, black) are drawn for p = lO 10-0 5 - 7 g/cm 3 with j — 1,2, ■ ■ ■ , 14, and the magnetic field lines (thin, 
green) are for A v = ^4i,9,max(l — 0.05j) with j = 1, 2, • • • , 19 where yl^.max denotes the maximum of A v . 



ner. First, we define a relativistic gravitational poten- 
tial $ c ee 1 - a c ($ c > 0), which is w 0.006 at t = 
for all the models chosen in this work. Since $ c is ap- 
proximately proportional to M/R, ^J 1 can be used as a 
measure of the characteristic length scale. From t = to 
the time at which $ c = 0.02, we set N = 420. The value 
of L is chosen so that the equatorial radius is initially 
covered by 400 grid points. At $ c = 0.02, the charac- 
teristic stellar radius becomes approximately one third 
of the initial value. At this time, the first regridding is 
performed. The grid spacing is cut in half and the grid 
number is increased to N — 700. Next, we set A" = 1180 
for 0.04 < $ c < 0.08, N = 1900 for 0.08 < $ c < 0.16, 
and N — 2500 for $ c > 0.16 (the minimum value of $ c is 
0.25 for all the models we consider). In this treatment, 
the total discarded fraction of the baryon rest mass is 
- 1% for model A, - 2% for model B, and - 2.5% 
for model C with (Ti, T 2 , T th , p nuc ) = (1.3, 2.5, 1.3, 
2 x 10 14 g/cm 3 ). In the following, we refer to this reso- 
lution as the standard resolution. 

To check the convergence of our numerical results, we 
also perform simulations using higher and lower grid res- 
olutions for model A2. In the higher-resolution case, 
the value of N is changed as follows: A~ = 500 for 
$ c < 0.02, N = 840 for 0.02 < $ c < 0.04, N = 1420 
for 0.04 < $ c < 0.08, N = 2268 for 0.08 < $ c < 0.16, 
and N = 2540 for <I> C > 0.16. In this case, the equa- 
torial radius is covered by 480 grid points initially. In 
the lower-resolution simulations, the equatorial radius is 
initially covered by 300 and 240 grid points and the val- 
ues of N are changed as follows: A^ = 320 and 260 for 
$ c < 0.02, N = 540 and 432 for 0.02 < $ c < 0.04, 
A" = 900 and 720 for 0.04 < $ c < 0.08, N = 1440 and 
1200 for 0.08 < $ c < 0.16, and N = 1900 and 1520 for 



$ c > 0.16. These resolutions are referred to as the high, 
middle, and low resolutions, respectively. 



We also performed several simulations using multiple 
transition fisheye coordinates [34] . This technique allows 
us to obtain high resolution in the central region with a 
relatively small grid number. The details of our fisheye 
implementation are discussed in Appendix A. We have 
confirmed that the results obtained using the fisheye co- 
ordinates agree with those obtained using the original 
coordinates. 



Simulations for each model with the standard resolu- 
tion are performed for about 150,000 time steps. This 
corresponds to the physical time of ~ 150 ms (~ 35 ms 
after core bounce). At this time, most of the interesting 
MHD processes have occurred and the system is settling 
down to a stationary state. For the code of [31], the re- 
quired CPU time for one model is about 20 CPU hours 
using 48 processors of FACOM VPP 5000 at the data 
processing center of National Astronomical Observatory 
of Japan, about 60 CPU hours using 8 processors of NEC 
SX8 at Yukawa Institute of Theoretical Physics (YITP) 
in Kyoto University, and about 120 CPU hours using 
8 processors of NEC SX6 at ISAS, JAXA. In order to 
check the validity of our results, several simulations were 
repeated using the code of [30]. We have confirmed that 
the results obtained with these two codes agree. Thus, 
in the following section, we mainly present the numerical 
results from the code of [31]. 
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VI. NUMERICAL RESULTS 



A. Dynamics 
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FIG. 2: Evolution of the central lapse and central density for 
models AO (dotted curves), Al (long-dashed curves), A2 (solid 
curves), A3 (dashed curves), and A4 (dot-dashed curves). 



As described in, e.g., [55, 56, 59], the collapse of a 
rotating stellar core for which the initial T/\ W\ is not too 
large (i.e., < 0.01) and the degree of differential rotation 
is small (A <; 0.5), can be divided into the following three 
phases. The first is the infall phase, in which the collapse 
sets in due to the sudden softening of the EOS. During 
this phase, the central density monotonically increases 
until it reaches nuclear density. The interior of the core 
that collapses nearly homologously is referred to as the 
inner core. The duration of the infall phase (i.e., the 
interval between the onset of the collapse and the time 
at which the central density reaches a maximum) is of 
order p c . 

The second phase is the bounce, which sets in when 
the densities in the central region exceed nuclear density 
p nuc - At this phase, the infall of the inner core decelerates 

— 1/2 

on a typical timescale of a few ms (~ 10p n uc )• Because 
of its large inertia and kinetic energy, the inner core does 
not settle down to a stationary state immediately but 
overshoots and rebounds, driving outgoing shocks at the 
outer edge of the inner core. 

The third phase is the ring-down. The bounce oc- 
curs when the central density reaches ~ 3p nuc due to 
a sudden stiffening of the EOS (if the centrifugal force 
is sufficiently small at the time that the density of the 
inner core exceeds nuclear density). In this case, the in- 
ner core oscillates quasi-radially for about 10 ms, and 
then settles down to a stationary state. In the outer re- 
gion, on the other hand, shock waves propagate outward, 
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FIG. 3: Evolution of the angular velocity of the resultant 
PNSs as a function of cylindrical radius in the equatorial plane 
for models AO, Al, A2, and A4. The long-dashed, dashed, 
dotted, and solid curves denote f2 at 118.4, 127.9, 136.7, and 
149.1 ms for model AO, at 118.5, 130.6, 140.3, and 150.6 ms for 
model Al, and at 118.3, 131.3, 137.8, and 150.8 ms for model 
A2, and at 118.6, 128.0, 137.2, and 149.5 ms for model A4. 
The dot-dashed line segments denote the slope of f2 oc -uj~ 3 ^ 2 . 



sweeping through infalling material from the outer enve- 
lope. With our simple treatment of the microphysics, the 
shock propagates all the way through to the surface of 
the outer core, which has a radius ~ 1000 km [69]. 

Even in the presence of magnetic fields, these general 
features are not significantly modified as long as the field 
is not extremely strong (i.e., as long as the magnetic en- 
ergy is smaller than the rotational kinetic energy). In 
Fig. 2, we show the evolution of the central lapse and 
central density for models A0-A4. The results for mod- 
els with non-zero magnetic fields (models A1-A4) are 
qualitatively very similar to those for model AO, which 
has no magnetic field. However, the magnetic effects lead 
to some quantitative changes. Figure 3 shows the evolu- 
tion of the angular velocity for the newly-formed PNSs 
as a function of the cylindrical radius in the equatorial 
plane for models AO, Al, A2, and A4. For model AO (see 
Fig. 3, upper panel), f2 relaxes to a stationary profile 
by ~ 10 ms after formation of the PNS. The remaining 
panels show the evolution of f2 in the presence of mag- 
netic fields (up to ~ 10 16 G). In these cases, the PNS 
slows down monotonically due to outward angular mo- 
mentum transport driven by magnetic braking, the MRI 
and the MHD outflow (see Sees. VI B and VIC for de- 
tails). In particular, for model A2, in which the growth 
of the magnetic field by winding saturates shortly after 
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FIG. 4: Evolution of the density profile of the PNSs as a 
function of the cylindrical radius in the equatorial plane for 
models AO (upper panel) and A2 (lower panel). The profiles 
are given at times corresponding to those of Fig. 3. The dot- 
dashed line segments denote the slope of p oc zu~ 4 . 



bounce (see Fig. 6), the central spin period doubles in the 
first ~ 20 ms following the bounce. For model Al, the 
magnetic field is amplified more gradually (see Fig. 6). 
In this case, does not decrease as rapidly as for model 
A2. For model A4, the PNS also spins down monoton- 
ically. This model has a different initial magnetic field 
profile from those of models A1-A3. This difference is 
reflected in the evolution of the angular velocity profile, 
which is not as smooth for A4 in the outer region of the 
PNS at late times. Since the magnetic field is stronger in 
the outer layers for model A4, the MRI is more vigorous 
there and the f2 profile is thus not as smooth (see also 
Sec. VI B). 

Figure 4 shows the evolution of the density profiles 
in the equatorial plane for models AO and A2. For AO, 
the density profile quickly relaxes to a stationary state, 
while it changes with time for A2 due to the angular mo- 
mentum transport. In accord with the spindown shown 
in Fig. 3, the central density gradually increases (see 
Fig. 2). Since material is ejected in an MHD-driven out- 
flow, the density in the surface region decreases. A more 
detailed discussion and a qualitative explanation of the 
main mechanisms driving angular momentum transport 
and the outflow are given in Sec. VIC. 

In Figs. 5(a) and (b), we show the various energies 
(E- lnt , Bheat) T mt , and Eem) as functions of time for 
models A0-A4, B, and C. All of the energies increase 
during the infall phase and reach their maxima at the 
bounce. Following the bounce, E- mt and E^eat relax to 
values of a quasistationary state which change very slowly 
(see Fig. 5(a)). On the other hand, T rot rapidly decreases 



in the presence of the magnetic field, reflecting the spin- 
down of PNSs by magnetic effects (see Fig. 5(b)). The 
rate of this decrease is larger for larger values of £em> as 
anticipated in Sec. HIE. 

Figure 5(c) shows the evolution ofE lnt , -Ehcat, Trot) and 
E~em for model A2 with different grid resolutions. With 
the middle and low resolutions, convergence is not well 
achieved. On the other hand, the results in the high and 
standard resolutions agree well, indicating that the stan- 
dard resolution is an appropriate choice for this problem. 

The simulations for models B and C are performed 
with the same initial magnetic field profile and field 
strengths as models A2 and Al, respectively (see Ta- 
ble II) , and the results for the various energies are shown 
together in Fig. 5. It is found that the evolutions for mod- 
els B and C proceed in the same qualitative manner as for 
models A2 and Al, respectively. Quantitatively, the ini- 
tial rotational kinetic energy is reflected in the evolution 
of T ro t and Eem- For model B, the maximum rotational 
kinetic energy is smaller than that for model A2. As a re- 
sult, the growth rate of the magnetic energy after bounce 
due to winding as well as the resulting maximum value 
reached are smaller. Thus, the magnetic energy achieved 
at saturation for a particular model depends on the rota- 
tional kinetic energy. This is expected since the magnetic 
energy is ultimately drawn from the energy stored in dif- 
ferential rotation. For model C, T ro t at bounce is larger 
than that of model Al because of the differential rota- 
tion at the onset of collapse and resulting larger rotation 
velocity of the PNS. As a result, the growth rate of -Eem 
and the maximum value reached are larger than those for 
model Al. 



B. Evolution of the magnetic field 

Figure 3 shows that for w ^ 10 km (i.e., in the central 
region of the PNSs), the angular velocity has a roughly 
flat profile. With so little differential rotation, the mag- 
netic field does not grow much, either by winding or the 
MRI. MHD effects have little influence on the evolution 
in this region. For w <; 10 km (i.e., near the surface 
and outside the PNSs), the angular velocity is approx- 
imately Keplerian (i.e., fl cx nj~ 3 / 2 ; see the dot-dashed 
line segments in Fig. 3). Because of the strong differen- 
tial rotation, the outer region with w ^ 10 km is subject 
to winding and the MRI (see Sees. IIIB and IIIC). In 
particular, the magnetic field is rapidly amplified in the 
region with vj ss 10-30 km, where both the angular ve- 
locity and the degree of differential rotation are large. 
We note that strong differential rotation results even for 
rigidly rotating progenitors. 

In Fig. 6, we show the evolution of the maximum val- 
ues of \B l \ for models Al, A2, and A4. Note that, 
in our setup, B x and B z are the poloidal field compo- 
nents (B x — _B ro ) and B y is the toroidal component 
(B y = wB v = B T ) since we assume axisymmetry and 
solve the equations in x-z plane. To demonstrate the ef- 
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FIG. 5: (a) Evolution of Bint and E^cat for models AO (dotted curves), Al (dashed curves), A2 (solid curves), A3 (dot-dashed 
curves), A4 (long-dashed curves), B (dot-long-dashed curves), and C (long-and-short dashed curves), (b) The same as (a) but 
for T ro t, and Eem- (c) Evolution of E- m t, -Ehcat, Trot, and Eem for model A2 with different grid resolutions. The dotted, solid, 
dashed, and dot-dashed curves denote the results in the high-, standard-, middle-, and low-resolution runs. 



fects of the grid resolution, the results for four different 
resolutions are displayed together for model A2. 

During the infall phase, the magnetic field is ampli- 
fied monotonically: The poloidal field grows due to the 
compression of matter, while the toroidal field grows pri- 
marily by winding (cf. Sec. Ill A). In the ideal MHD 
approximation, the magnetic field lines are frozen into 
the plasma, and hence, the poloidal field strength should 
grow approximately as p 2 / 3 (see Sec. Ill A). As shown 
in Fig. 6(c), this relationship holds for the maximum 
poloidal field components over several decades of growth 



in the density. Since B y grows by winding, its growth 
accelerates rapidly just before the bounce due to the 
rapid growth of \B X \. As shown in Fig. 6, a significant 
amplification indeed occurs for about 2-3 ms before the 
bounce at t 115 ms. Figure 6(c) also shows that the 
maximum value of the toroidal magnetic field increases 
slightly faster than p 2 / 3 during the infall phase. This 
qualitatively agrees with the estimate in Sec. Ill A. 

To provide an overview of the evolution, we display 
snapshots of density contours, velocity vectors, poloidal 
magnetic field lines, and magnetic pressure contours for 
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FIG. 6: (a) Evolution of \B l \ for model A2. The solid, dashed, long-dashed, and dot-dashed curves denote the results from 
the standard, high, middle, and low resolution runs for model A2. The enlarged panel for the evolution of |-B x | roa3t illustrates 
that the growth of the field strength by the MRI is well resolved only in the high and standard runs, (b) The same as (a) 
but for models Al (dotted curves), A2 (solid curves), and A4 (dashed curves) with the standard resolution, (c) l-B^lmax as a 
function of p m&K during infall phase for models Al (dotted curves), A2 (solid curves), and A4 (dashed curves) with the standard 
resolution. The dot-dashed line segments above the three curves denote the relation \B l 
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model A2 at selected time slices in Fig. 7. In Fig. 8, 
the evolution of the magnetic field strength as a function 
of the cylindrical radius is shown for selected coordinate 
values of z. (We do not plot B x and B y on the equa- 
torial plane since they vanish there by symmetry.) We 
find that, after the bounce at t ~ 115 ms, shock waves 
propagate outward. Material behind the shocks moves 
outward in an anisotropic manner due to the rotation 
and to the anisotropic density profile of the collapsing 



star. This matter outflow is strongest along the rotation 
axis. 

For t <; 120 ms, differential rotation modifies the pro- 
file of the magnetic field and amplifies the magnetic pres- 
sure (see Figs. 6, 7, and Fig. 8). As Fig. 6 shows, the 
toroidal field continues to grow after bounce because of 
the differential rotation near the surface and outside of 
the PNS. The magnetic pressure monotonically increases 
for 120 ms & t <, 130 ms in the region w = 10-30 km 
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FIG. 7: Density contour curves and velocity vectors for model A2 (first and third rows). The contours are drawn for p — 
10 i5-o.5i g / cm 3 (j _ i_i2). (The red, blue, magenta, cyan, and yellow curves denote p = 10 14 , I0 13 , I0 12 , 10 11 , and I0 10 g/cm 3 , 
respectively.) The scale of the velocity is shown in the upper-right corner. The second and fourth rows show the poloidal 
magnetic field lines (green) and contours of the magnetic pressure (thick red, blue, and normal black curves) at corresponding 
times. The poloidal magnetic field lines are drawn as contours of A v , with levels given by A v = (1 — Q.li)A v , max (i = 0-9). 
Here, v4 v , max is the maximum value of A v at each time slice. The contour curves of the magnetic pressure are drawn for 10 30 
(very thick red curves), 10 29 (thick blue curves), and 10 28 dyn/cm 2 (thick black dotted curves). 



and z <J 10 km. This is reflected in the third through 
fifth snapshots in Fig. 7, which show that the region with 
fmag > 10 30 dyn/cm 2 expands. We note that the am- 
plification at z — is very small since we impose the 
symmetry condition that B x and B v must vanish on the 
equatorial plane. Comparison of Figs. 5 and 6 illustrates 



that the increase of the magnetic energy is primarily due 
to the growth of the toroidal field. This growth stops 
when the magnetic energy reaches about 10-20% of the 
rotational kinetic energy irrespective of the initial field 
strength and profile (cf. Fig. 5). Figure 6 shows that the 
toroidal field growth can be followed well even for the 
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lower grid resolutions. 

The growth of magnetic energy saturates when back- 
reaction from magnetic braking becomes significant. This 
occurs at a time comparable to the Alfven crossing 
time [13] t sa t ~ t\, where tA is given by Eq. (52). For 
the models considered here, t\ 30 ms. From the nu- 
merical results for model A2, we see that £ sa t ~ 10 ms 
(see Fig. 6). 

The toroidal magnetic field would not achieve the sat- 
uration strength described above if significant magnetic 
reconnection were to operate before magnetic braking 
could take effect. However, the velocity of fluid enter- 
ing the dissipation layer is expected to be a small frac- 
tion (probably 10 _1 or 10~ 2 ) of the Alfven velocity (see, 
e.g. [74-76]). Simple estimates show that this expecta- 
tion is borne out in the solar magnetosphere and the 
geomagnetic tail of the earth [76]. Thus, the reconnec- 
tion timescale should, in general, be of order 10 to 100 
times longer than the Alfven timescale which governs the 
growth of the toroidal field. On the timescales of interest 
for PNS evolution considered here, reconnection will thus 
not play an important role. 

Following saturation, the toroidal field strength begins 
to decrease in the central region due to magnetic braking. 
This is reflected in the decrease of the magnetic pressure 
for vj ~ 10 km and z ~ 10 km seen in the sixth and 
seventh snapshots of Fig. 7. Magnetic braking also leads 
to outward transport of angular momentum by Alfven 
waves. This is clearly seen in the middle panel of Fig. 8; 
\B V \ decreases for w <, 15 km but increases for tu ^ 15 
km. The middle and lower panels of Fig. 3 also show 
that the angular velocity in the central region decreases 
while it increases for 10 km <, w ^ 50 km for models 
Al and A2. This angular momentum transport drives 
matter outflow at relatively low latitude (see for t ^> 130 
ms in Fig. 7). However this outflow is not as strong as 
that along the rotation axis [77] . 

At bounce, the compression stops and so does the rapid 
global growth of the poloidal field. However, the poloidal 
field in the outer region of the PNSs (see Fig. 7) continues 
to grow after the bounce by the MRI. Figure 7 indicates 
that the poloidal field lines for vd sa 10-30 km and z — 
10-50 km are highly distorted for t ^> 130 ms. Evidence 
of the MRI is also seen in the first and third panels of 
Fig. 8 which show that initially smooth poloidal fields 
become highly inhomogeneous and the amplification is 
locally enhanced [77] (the local amplification can be seen 
clearly in the region R ~ 35 km in the B x plot in Fig. 8). 

The growth l-B^lmax shown in Fig. 6(a) demonstrates 
the importance of resolution in capturing the MRI. In 
this figure, we plot results for four different resolutions. 
In the low and middle resolutions, the maximum value of 
\B X \ increases until the bounce and holds steady there- 
after. On the other hand, there is significant growth 
after the bounce in the standard and high resolution 
runs. Note that the wavelength of the fastest-growing 
mode of the MRI, A max , is a few km for B ~ 10 14 G, 
p ~ 10 13 g/cm 3 , and fi ~ 10 3 rad s" 1 [cf. Eq. (60)] 
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FIG. 8: Evolution of \B X \ and \B V \ atzRi 14.2 km and \B Z \ 
in the equatorial plane as a function of the cylindrical radius 
for model A2. The solid, dotted, and dashed curves denote 
the data for t = 150.8, 137.8, and 118.3 ms. The thick solid 
curve for \B Z \ is the result at z fa 10 km at t = 150.8 ms. 



which are typical values in the outer region of the PNS. 
In the standard resolution, the grid spacing is ~ 0.1A max , 
and the MRI is marginally resolved, but it is unresolved 
for the low and medium resolutions. This illustrates that, 
in MHD simulations, integrations with several grid reso- 
lutions are essential for correctly identifying the physical 
mechanisms at work. 

Figures 6(a) and (b) show that |-B 2 | max does not in- 
crease significantly. This does not imply that |-B 2 | does 
not change due to the MRI. Indeed, we find that the z- 
component is amplified in the outer regions of the PNSs 
(see Fig. 8). However, these local fluctuations are not as 
large as B z in the central part of the PNSs (where the 
z-component dominates and is fairly constant following 
the PNS formation). 

The amplification of the magnetic field by the MRI is 
seen at ~ 10 ms after the bounce. The estimated growth 
time scale of the MRI is t M m » 4/3fi [cf. Eq. (59)] 
which is a few milliseconds for w — 10 30 km. Our 
simulation indicates that the actual iMRi is somewhat 
longer than this estimate. The tendency for the linear 
analysis to underestimate the MRI growth time is also 
seen in our previous results for the evolution of mag- 
netized neutron stars [33]. The linear analysis may be 
inaccurate by a factor of several in the present context 
since the MRI is treated as a purely local phenomenon, 
assuming a uniform background state over length scales 
much longer than the wavelengths of the perturbations. 
However, since the expected value of A max is only one 
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order of magnitude smaller than the characteristic ra- 
dius, this assumption may not hold. To understand the 
growth timescale of the MRI in the PNS, a global anal- 
ysis is necessary. Indeed, a global stability analysis for 
magnetized accretion disks shows that the local analysis 
overestimates the growth rate [78]. 

As in the case of the toroidal field, the growth of the 
poloidal field saturates when l-B^lma* reaches ~ 2-3 x 10 16 
G. The saturation may be partly due to the limited grid 
resolution (i.e., due to numerical diffusion), but this effect 
does not seem to be very large since the value of l-B^max 
is approximately the same for the high and standard res- 
olutions (see Fig. 6). 

The contribution to the growth of the magnetic field 
strength by winding is significantly larger than the con- 
tribution from the MRI. Indeed, we find that the aver- 
age value of the poloidal field is smaller than that of the 
toroidal field. Because MRI is a local instability, it does 
not increase the global magnetic energy as significantly 
as winding. 

After the saturation of the field growth by winding and 
the MRI, a collimated and nearly stationary magnetic 
field with \B Z \ ^> \B rxl \ is formed near the rotation axis 
(see the last snapshot of Fig. 7). For the high latitude 
regions with z <; 50 km and with w ^ 10 km, the abso- 
lute value of the toroidal component \B V \ is as large as 
\B Z \, indicating a helical magnetic field. (On the other 
hand, very near the rotation axis, \B Z \ is much larger 
than The field strength at t ~ 150 ms is £ 10 15 G 

at z = 100 km. Since the degree of differential rotation 
is small near the rotation axis, this magnetic field config- 
uration is stable. A nearly stationary outflow is driven 
along the collimated field lines. This is probably due to 
the magneto-centrifugal mechanism (cf. Sec. HIE), and 
we discuss this possibility in more detail in Sec. VI C. In 
contrast to the region near the axis, the toroidal field is 
much stronger than the poloidal field at large cylindrical 
radius. 

We find that the evolution path of the magnetic field 
depends quantitatively (but not qualitatively) on the ini- 
tial field strength and profile. In the right panel of Fig. 6, 
we show the evolution of each component of the magnetic 
field for models Al, A2, and A4 (see also Fig. 5). For Al, 
the initial magnetic field profile is the same as that for 
A2 while the initial strength is weaker. The evolution of 
the magnetic field is similar for these two models except 
that, for Al, it takes longer to amplify the magnetic field 
because of the smaller initial value of £? ro . 

Model A4 has a different initial magnetic field profile 
from that of models Al and A2. Since the magnetic field 
is initially distributed over a larger range of cylindrical 
radii for A4 (see Fig. 1), field growth by compression is 
more efficient for A4. Although the maximum values of 
the magnetic field are initially identical for models Al 
and A4 (see Table II), those at bounce are larger for 
model A4. Model A4 also has a different field configura- 
tion after bounce than Al and A2. In Fig. 9, we show 
snapshots of the density contours, the velocity vectors, 



the magnetic field lines, and contour curves of the mag- 
netic pressure after the formation of the PNS for model 
A4. By comparison between Fig. 7 and Fig. 9 at t w 116 
ms, it is found that the fraction of the z-component of 
the magnetic field is larger for model A4 than for model 
A2. This strong vertical field is advantageous for induc- 
ing the MRI. In addition, a stationary collimated field is 
formed more quickly for A4 than for A2. This results in 
a stronger MHD outflow (see Sec. VIC). 



C. MHD outflow 

Both winding and the MRI increase the magnetic stress 
in the outer regions of PNSs which have large degrees of 
differential rotation. This stress induces an MHD out- 
flow, particularly near the rotation axis (see velocity vec- 
tors of Figs. 7 and 9 for t <; 130 ms). To prove that 
this is indeed due to the magnetic stress, in Fig. 10, we 
show evolution of the density contour curves and veloc- 
ity vectors for model AO, which has no magnetic field. 
In this model, the PNS relaxes to a stationary state for 
t <; 130 ms. Comparing Figs. 7, 9, and 10, it is ev- 
ident that the MHD outflow is driven by magnetic ef- 
fects. One also sees that the density profile changes due 
to the matter outflow (especially near the rotation axis) 
and that the oblateness of the PNSs for models A2 and 
A4 is smaller than for model AO at t w 150 ms due to 
the angular momentum loss (compare the density con- 
tour curves of p — 10 12 g/cm 3 for the three models in 
the final snapshots). 

We posit three possible sources for the MHD outflow: 
the magneto-spring mechanism, turbulence induced by 
the MRI, and the magneto-centrifugal mechanism. In 
the early phase, winding amplifies the toroidal field and 
the magneto-spring mechanism works efficiently to drive 
the MHD outflow. This outflow blows off matter pref- 
erentially in the z-direction. However, it is not sharply 
collimated; the half-opening angle is ~ 45 degrees. Tur- 
bulent motion caused by the MRI also induces an outflow. 
In contrast to the magneto-spring mechanism, however, 
the MRI-driven outflow is incoherent. 

After the toroidal magnetic field saturates, a station- 
ary collimated, helical magnetic field forms around the 
rotation axis. Matter then flows out along the collimated 
field lines for t ^ 150 ms in case A2 and for t <; 140 ms in 
case A4, as shown in Figs. 7 and 9. This suggests that the 
magneto-centrifugal mechanism (see Sec. HIE) is oper- 
ating. This interpretation can be inferred indirectly from 
two facts: (i) the MHD outflow is continuously ejected 
from the PNS, even after the growth of the toroidal field 
by winding saturates, and (ii) the outflow is absent in 
the region very close to the rotation axis. This indicates 
that the outflow is driven along the collimated magnetic 
field with an appropriate inclination angle between the 
rotation axis and the poloidal field lines, as is necessary 
for magneto-centrifugal launching [66, 67]. 

Further evidence is found by confirming that the quan- 




FIG. 9: Evolution of the density contour curves and velocity vectors for model A4. The contours and velocity vectors are drawn 
in the same manner as in Fig. 7. 
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FIG. 10: Evolution of the density contours and velocity vectors for model AO. The contours and velocity vectors are drawn in 
the same manner as in Fig. 7. 



143 ms 



t- - - 1 






- ^ 












7 




-1 -, \ 



tities k = p*v p /B p , which should be constant along 
the field lines in a stationary state are actually con- 
stant (see [66, 79] and Appendix B). Here v p and B p 
denote the poloidal components of v l and B l . Note that 
the definition of k accounts for general relativistic cor- 
rections through p* and B p . In Fig. 11, we show the 
contour curves for k together with the magnetic field 
lines at t = 152.5 ms for model A4. It is found that 
the two sets of curves overlap in the region of the colli- 
mated magnetic field lines, verifying that the structure 
of the magnetic field is almost stationary in this region. 
On the other hand, the two sets of curves do not overlap 
for the turbulent region far from the rotation axis. 

To give a better overall picture of the outflow, we show 
in Fig. 12 the density contours, velocity vectors, and 
poloidal magnetic field lines at t = 152.5 ms for model 



A4. In this figure, a larger region than that shown in 
the last panel of Fig. 9 is displayed. For r ^ 200 km, a 
weakly-anisotropic outflow is found. However, the veloc- 
ity of the matter near the rotation axis is relatively large. 
For r 200 km, on the other hand, the strong outflow 
is seen only along the collimated magnetic field in the 
vicinity of rotation axis. The greatly enhanced strength 
of the outflow along the collimated field lines points to 
magneto-centrifugal launching as a likely driver of the 
outflow. Our results suggest that jets may be launched 
during supernova explosions if the collimated magnetic 
fields found here are generic [6, 7, 28]. For r <S 200 km, 
a slower, incoherent flow pattern is also seen, reflecting 
the continuous driving of irregular matter motions by the 
MRI. 

As mentioned in Sec. HIE, the magneto-centrifugal 
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FIG. 11: Contour of k = p*v p /B p (solid curves) and mag- 
netic field lines (dotted curves) at t = 152.5 ms for model 
A4. The magnetic field lines are drawn in the same man- 
ner as in Fig. 7. The contours of k are drawn for k — 
6 x I0 4 - 4i g/(cm 2 sG)(t = 0-4). We note that these con- 
tours are drawn only for the values of k found for material on 
the collimated field lines. 



launching mechanism requires a mechanism for ejecting 
material outward at the base. This is particularly im- 
portant if the opening angle of the collimated field line 
is small [66, 67]. In the present case, the magnetic field 
is significantly perturbed in the vicinity of the PNSs by 
the MRI. Turbulence and/or irregular Poynting flux gen- 
erated by the MRI could serve to inject material into the 
magneto-centrifugal wind. 

In addition to rest mass, energy and angular momen- 
tum are also carried away from the PNS by the outflows. 
In Fig. 13 we show Fm, F e , and Fj evaluated at r = L/4 
(w 224 km for models A0-A4, w 143 km for model B, 
and w 161 km for model C) as a function of time for 
models A0-A4, B, and C. We also evaluated the fluxes 
at other radii and found that, aside from a time shift, the 
results depend weakly on the radius. 

In the early ringdown phase in which shocks propagate 
outward, the loss of mass, energy, and angular momen- 
tum from the system is due primarily to the matter out- 
flow (as opposed to the Poynting outflow) irrespective of 
the presence of magnetic field. Hence Fm, F e , and Fj 
are roughly independent of the field strength (compare 
results for model AO, Al, A2, and A4). In the absence of 
the magnetic field (model AO; dot-dashed curves), how- 
ever, the outflow rates decay in an exponential manner, 
leaving a stationary PNS. On the other hand, in the pres- 
ence of the magnetic field, the fluxes remain nonzero, and 
the PNSs continuously lose mass, energy, and angular 
momentum. The evolution of the fluxes for model Al 
show that the dominance of the MHD-driven outflow 



over the shock-driven outflow is delayed for ~ 10 ms 
after bounce, i.e. after the toroidal field becomes suffi- 
ciently strong. (Note that for model Al, the growth of 
the toroidal field by winding saturates at the relatively 
late time of ~ 30 ms after bounce. This is consistent with 
the fact that the MHD-driven outflow becomes dominant 
later for this model than for A2-A4.) 

The strength of the outflow depends on the initial mag- 
netic field strength (compare the results of models Al 
and A2). This is due to the difference in the resulting 
ratio of the poloidal field strength to the toroidal field 
strength C B for the PNS [see Eq. (64)]. The toroidal 
fields are wound up to similar saturation strengths in 
the various models regardless of the initial poloidal field 
strength (see Fig. 6). On the other hand, after bounce, 
the poloidal field is not uniformly amplified, though the 
MRI increases the strength locally. As a result, the value 
of Cb is smaller for model Al and consequently so are 
the values of Fm, F e , and Fj. 

In addition to the magnitude of the initial seed mag- 
netic field, the strength of the outflow (and particularly of 
the magneto-centrifugal outflow), also depends strongly 
on the initial field profile. The relative size of the z- 
component of the magnetic field is larger for model A4 
than for model A2 (compare Figs. 7 and 9). The larger 
z-component is favorable for inducing the MRI and the 
resulting MHD outflow, although the magnetic energy of 
the PNSs for models A2 and A4 are nearly identical (see 
Fig. 5). 

Although the outflow rates depend on the initial mag- 
netic field profile and strength, the ratio of Fj / Fm ~ 1.5- 
2.5J/M is always larger than J/M for all cases consid- 
ered here. Thus, the MHD outflow carries away material 
with a larger specific angular momentum than that of 
the PNSs, leading to spindown. Figure 13 shows that 
the spindown time scale Jq/Fj (Jo is the initial value 
of J) is ~100-300 ms for models Al, A2, and A4, and 
is approximately proportional to Cg 1 [see Eq. (64) and 
note that the value of Cb for A2 is about twice of that 
for Al]. Since Cb <; 0.1 is indicated by our simulations, 
a large amount of the angular momentum of the PNSs 
is carried away within ~ 1 s after the onset of the MHD 
outflow. 

In Fig. 14, we show the ratios of Fe. em /Fe and 
Fj,em/Fj as functions of time for models Al, A2, A4, 
B, and C. The plot indicates that the outflow is largely 
Poynting dominated. The ratio of the EM flux to the 
total flux clearly depends on the initial magnetic field 
strength and profile. Comparing the results shown in 
Fig. 14, F E . EM /F e and Fj. em /Fj are larger for larger 
values of Fm and Fj [80]. 

The large relative fraction of the Poynting flux at radii 
around ~ 200 km also implies that, during the outward 
propagation of the magnetic energy, conversion to mat- 
ter kinetic energy is suppressed in our model. In realistic 
supernovae, shock waves formed at bounce do not prop- 
agate outward unimpeded but rather stall at a radius 
~ 150-200 km, forming a standing shock [70-73]. Thus, 
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FIG. 12: The density contours and velocity vectors (left) and the poloidal magnetic field lines and velocity vectors (right) at 
t — 152.5 ms for model A4. The density contours are drawn for p = 10 15_! g/cm 3 (i — 1-7). The velocity vectors and poloidal 
magnetic field lines are drawn as in Fig. 7. 
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FIG. 13: F M , F e , and Fj at r « 224 km for models AO (dot- 
long-dashed curves), Al (dashed curves), A2 (solid curves), 
A4 (long-dashed curves), B (dotted curves), and C (dot- 
dashed curves). The abbreviation 'foe' denotes 10 51 ergs, 
while Jo is the initial value of J. 
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FIG. 14: F E , em /F e and Fj,em/Fj for models Al (dashed 
curves), A2 (solid curves), A4 (long-dashed curves), B (dotted 
curves), and C (dot-dashed curves). 



the MHD outflow driven from the PNS will eventually 
hit the stalled shock, where the Poynting energy flux as 
well as the kinetic energy of the matter outflow may be 
converted to thermal energy. The total energy from the 
outflow during the 10 ms period following the relaxation 
of Fm to a roughly constant value is <; 10 50 ergs for all 
of the models. If this energy were converted to thermal 



energy at the stalled shock, it could play a significant 
role in driving the supernova explosion, as pointed out in 
[9, 27, 28]. 

In Fig. 15, the ratio of Fj to -Eem is shown for mod- 
els Al, A2, A4, B, and C. By Eq. (66), this ratio should 
be of order Cb if the MHD outflow theory [66] holds in 
the vicinity of the PNS. Figure 15 shows that, in the early 
phase (t <i 130 ms) in which matter outflow driven by 
shocks is dominant, the value Fj / -Eem decreases steeply. 
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FIG. 15: Fj/E E m for models Al (dashed curve), A2 (solid 
curve), A4 (long-dashed curve), B (dotted curve), C (dot- 
dashed curve). 



However, once the MHD outflow dominates, the ratio re- 
laxes to an approximately constant value of ~ 0.1-0.2 
for all of the models. (Our simulations indicate that 
C B = \B P \/\B T \ ~ 0.1-1 for these models.) We also find 
that its value for model Al is about half as large as the 
values for models A2 and A4, reflecting the smaller value 
of Cb for model Al. Thus, our results agree qualita- 
tively with the MHD outflow theory. The MHD outflows 
found here are ultimately powered by rotation. Thus, the 
outflow is expected to weaken as the PNSs spin down, al- 
though at least in the first few tens of milliseconds after 
the saturation of the field growth, the strength is approx- 
imately constant as shown in Fig. 15. 



D. Outflow-induced spindown of the PNSs 

Assuming that the PNSs spin down due to the MHD 
outflow, the loss rate of the angular momentum is de- 
scribed by Eq. (66) with 77* = O(l). The spin angu- 
lar momentum of the PNS is approximately written by 
</* = 7*r2* where P is a moment of inertia and fL is an 
average angular velocity. After the saturation of the mag- 
netic field growth, Pem is approximately equal to (Trot, 
where C is around 0.1-0.2 according to our results (see 
Fig. 5), and where T rot is approximately I*Q%/2. Gath- 
ering these results, we have a relation for the evolution 
of the angular momentum of the PNS: 



d ( I* f2 ^ 
dT 



Assuming that P and Cb are constant, we then have 

n^C B , (75) 



PP. 
~df 



d /2vr 
dt \ fL 



where P* = 27r/fL- Equation (75) suggests that the 
rotational period increases linearly with time as 

P* » P +7rCry*Cflt 

where Po is the initial value of P*. The value of Cb 
during the stationary MHD outflow phase is somewhat 
unclear for a real system. Our present result, however, 
indicates that \B P \ ~ \B T \ near the rotation axis, giving 
C B ~ 0.1-1. 
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FIG. 16: P* as a function of time for model A4. 



In Fig. 16, we show the evolution of P* for model A4, 
in which the PNS and magnetic field lines reach a quasis- 
tationary state for t > 140 ms. To obtain P», we compute 
fL from 



p*nd 3 x 



p>10" 3 p„ 



p*d 3 x 



p>10" 3 p„ 



.(77) 



I^VI^Cb- 



(74) 



The figure shows that the average spin period increases 
approximately linearly in time for t J> 140 ms. This 
agrees with Eq. (76). For this case, we find dP^/dt sw 

0. 013. Thus the order of magnitude for the spindown 
rate is in agreement with Eq. (75). 

The rapid spindown will continue until the rotation 
in the vicinity of the PNS becomes sufficiently uniform, 

1. e., until the magnetic field relaxes approximately to a 
stationary state. If the magnetic field is amplified until 
saturation with £ ~ 0.1 shortly after bounce, and if the 
degree of differential rotation remains high for ~ 1000 s 
after the saturation, the rotation will slow to a period 
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longer than ~ 1 s. For poloidal magnetic field strengths > 
10 15 G, the resulting neutron star will be a magnetar [3, 
8]. However, we note that these results may be modified 
when the PNS cools, which occurs on a timescale (~ 20 s) 
which is much longer than that of our simulations. 

We know that a class of young pulsars (Crab-like pul- 
sars) has rotational periods of several 10 ms [81]. The 
birth of such rapidly rotating young pulsars thus seems 
to require that the differential rotation in the vicinity of 
the PNS is quickly lost in order to avoid winding and the 
MRI, which lead to outflows and spindown. 

E. Gravitational waveforms 

In Fig. 17, we show gravitational waveforms for mod- 
els A0-A4. Gravitational waveforms are computed using 
the quadrupole formula described in Sec. II D. According 
to the classification system of [55, 56], the waveforms for 
these models are of type I, which is the most common 
type of waveform. The properties of type I waveforms 
may be summarized as follows: During the infall phase, a 
gravitational wave precursor is emitted due to the chang- 
ing quadrupole moment as the collapsing core flattens. 
In the bounce phase, spiky burst waves are emitted for a 
short time ~ 1 ms, and the amplitude and the frequency 
of gravitational waves reaches a maximum. In the ring- 
down phase, gravitational waves associated with several 
oscillation modes of the PNS, for which the frequencies 
are ~ 1 kHz, are emitted with amplitudes which are grad- 
ually damped due to shock dissipation at the outer edge 
of the PNS. 
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FIG. 17: A2 as a function of time for models AO (dotted 
curve), Al (dashed curve), A2 (solid curve), A3 (dot-dashed 
curve), and A4 (long-dashed curve). 

We find that the MHD outflow modifies the gravita- 
tional waveforms in the ringdown phase; the wave am- 



plitude gradually increases with time during this phase, 
although the waveforms are otherwise unchanged. That 
this secular increase is caused by the MHD outflow may 
be understood through the following simple analysis. 
Consider a matter outflow of mass m and velocity v in the 
z direction, and assume that the velocity changes slowly. 
The contribution to A2 from this outflow comes from I zz , 
for which the correction is 8A2 « 2mv 2 . The MHD out- 
flow is continuously ejected from the vicinity of the PNS, 
and hence, the total mass of the outflow increases with 
an approximately constant rate. We find that A 2 indeed 
increases roughly linearly with time. Furthermore, the 
magnitude is approximately given by 

SA 2 ^3(^-)(^-) 2 m, (78) 

which is in approximate agreement with the numerical 
results for the mass increase rate dm/dt = O.l-5M /s. 

The secular increase of the wave amplitude for model 
A4 is faster than for model A2 although the magnetic 
energies of the PNSs are nearly identical. This is because 
the MHD outflow is more efficiently driven for model A4 
as shown in Sec. VI C. 

As discussed in Sec. VIC, the MHD outflow from a 
newly-formed PNS will likely hit a standing shock at a 
radius ~ 150-200 km [70-73]. The value of m (and hence 
the linear growth of the amplitude) will then saturate 
in a few 10 ms. Nevertheless, the amplitude A2 could 
eventually reach a few meters, which is comparable to 
the amplitude of the spiky waves emitted at the bounce. 
The outflow signal may thus be detectable for an event 
within our Galaxy by the ground-based laser interfero- 
metric detectors, since the expected amplitude is ~ 10~ 20 
[see Eq. (28)]. Detecting gravitational waves in the first 
few tens of ms of the supernova collapse event may thus 
provide information about the anisotropic outflow. 



F. Comparison of codes 

To demonstrate agreement between the codes of [30] 
and [31], we plot some representative quantities in 
Figs. 18-20. Figure 18 compares the central density evo- 
lution and late-time density profiles for model A2 from 
both codes. The plot of central density shows agreement 
on the maximum density at bounce, the slow increase in 
the central density caused by the spindown of the PNS, 
and the post-bounce oscillations (see the inset of Fig. 18). 
(Results from the code of [30] were shifted forward in time 
by 0.35 ms to align the core bounce. This small differ- 
ence in the coordinate time likely results from the slightly 
different lapse prescriptions used in the two simulations.) 
The two codes also give the same principal features of the 
density profile at late times (e.g. t — 150.8 ms), including 
the falloff behavior for the density outside the PNS. The 
behavior of the magnetic field components, as compared 
in Fig. 19, also shows good agreement. Some differences 
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arise after the PNS forms due to the fact that the mag- 
netic field evolution is dominated by the turbulent effects 
of the MM. In this regime, the precise value of the max- 
imum magnetic field component depends on the details 
of the turbulent motions; small differences in local quan- 
tities can be amplified by turbulence. Finally we show 
the evolution of the individual components of the energy 
for model A2 (compare to Fig. 5), which also demon- 
strates good agreement between the two codes. We note 
that both the central density and magnetic field compo- 
nents grow by a factor of order 10 4 during the evolution 
of model A2. Given this large dynamic range, the agree- 
ment shown here is rather noteworthy. 
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FIG. 18: Results for model A2 obtained with the codes of [31] 
(solid lines) and of [30] (dotted lines), (top panel) Evolution 
of the central density for model A2. The inset shows the 
bounce and immediate post-bounce behavior, (lower panel) 
Comparison of the density profiles at t — 150.8 ms. The dot- 
dashed line segment shows a slope of p oc zxj~ 4 . The results 
from the code of [31] correspond to the standard resolution. 
For the results from the code of [30], a fisheye grid was used 
for the PNS evolution phase (see Appendix A for details). 



VII. SUMMARY AND DISCUSSION 

We have presented general relativistic numerical sim- 
ulations of magnetized stellar core collapse to a PNS in 
axisymmetry and have followed the subsequent evolution 
of the PNS by magnetorotational effects. The following 
is a summary of the results: 

1. The magnetic field is amplified during the collapse 
and the subsequent evolution of the PNS. During the in- 
fall phase, the poloidal field grows due to the compression 
of the matter in which the field is frozen. As the poloidal 
field strength grows, the amplification of the toroidal field 
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FIG. 19: Evolution of the maximum values of \B l \ for model 
A2. As in Fig. 18, results from the codes of [31] (solid lines) 
and of [30] (dotted lines) are compared. 
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FIG. 20: Individual components of the energy for model A2, 
defined as in Fig. 5. Results from the codes of [31] and [30] 
are shown as solid and dotted lines, respectively. 



by winding is accelerated. After bounce, the compression 
stops, and so does the rapid growth of the poloidal field. 
However, the differential rotation in the vicinity of the 
PNS induces further growth of the toroidal field through 
winding. The winding continues until the magnetic field 
drains away the kinetic energy stored in differential ro- 
tation (which is roughly 10-20% of the total rotational 
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kinetic energy). The ratio of the magnetic energy to the 
rotational kinetic energy thus saturates at approximately 
the same value regardless of the initial magnetic field pro- 
file and strength. The poloidal magnetic field in the outer 
region of the PNS also grows via the MRI. In contrast to 
magnetic winding, the MRI amplifies the field strength in 
a turbulent manner. The magnetic stress increases until 
it becomes ~ 10% of the rotational kinetic energy density. 
After the magnetic field saturates, we find a stationary, 
collimated helical magnetic field near the rotation axis. 

2. Because of the enhanced magnetic stress, an MHD 
outflow is launched from the vicinity of the PNS, par- 
ticularly near the rotation axis. Our results suggest that 
three mechanisms drive the outflow. One is the magneto- 
spring mechanism in which magnetic stress primarily due 
to the toroidal field blows off matter along the rotation 
axis (but without strong collimation), leading to the so- 
called 'tower-like' field structure. The second mechanism 
is MRI-induced turbulence, which leads to a weaker and 
less coherent matter outflow than the magneto-spring 
effect. The third mechanism is the Blandford-Paync 
magneto-centrifugal mechanism, which plays a dominant 
role after the saturation of the magnetic field growth and 
after the formation of a collimated helical magnetic field 
near the rotation axis. In the present context, matter 
is injected into the outflow by MRI turbulence and then 
flung out along the helical magnetic field lines. By these 
three mechanisms, the energy carried away by the MHD 
outflow in the first 10 ms after the saturation of the 
toroidal magnetic field growth may be i> 10 50 ergs for 
a PNS of rotation period ~ 1-2 ms and for a typical ra- 
tio of the poloidal to toroidal field strength of Cb <; 0.1. 
3: The MHD outflow carries away material with large 
specific angular momentum, and hence, plays a crucial 
role in spinning down the PNSs. We find that the angular 
momentum loss rate is ~ 0.1-0.2_EemC_b as long as differ- 
ential rotation persists in the vicinity of the PNS. Since 
£em is ~ 10-20% of the rotational kinetic energy T TOt 
after the field saturates, we have \J\ ~ 0.01T rot CB- This 
implies that the angular momentum (and rotational ki- 
netic energy) may decrease rapidly as show in Sec. VI D. 
If the differential rotation remains strong enough in the 
first <~ 1000 s after the magnetic field saturates, the spin 
period of the PNS increases to 1 s. 

In turbulent flows resulting from the MRI, small eddies 
would be formed, dissipating kinetic energy into ther- 
mal energy and transporting angular momentum out- 
ward. The thermal energy generated in this process has 
been suggested as a power source for a supernova explo- 
sion [9] . However, we find that the role of the turbulence 
is likely less important as a source of energy than the 
MHD outflow. This may be partly due to our choice of 
the initial seed magnetic field. In the present work, the 
poloidal magnetic field at the time of PNS formation is 
high enough that the toroidal field is rapidly amplified 
and drives a strong MHD outflow in a few tens of ms af- 
ter the bounce. For a weaker initial poloidal field, more 
time would be required to amplify the toroidal field [see 



Eq. (51)] and, hence, to produce a strong MHD out- 
flow. In contrast, the MRI grows on the same timescale 
(<~ lOil -1 ) even for very weak fields. Since the wave- 
length of the fastest-growing MRI mode scales with the 
field strength, very weak fields will result in turbulence 
composed of small eddies, for which the typical size is 
much smaller than the stellar radius. This small scale 
turbulence converts kinetic energy into thermal energy, 
as modeled in [9]. In this case, the magneto-spring and 
magneto-centrifugal mechanisms may not play a signifi- 
cant role and the collimated magnetic field may not form. 
Unfortunately, simulations for such weak initial fields are 
currently unable to capture this behavior, since the wave- 
length of the MRI is much too short to be resolved. To 
determine whether such processes indeed occur, a much 
more powerful supercomputer will be required. 

We do not observe the magnetic field amplification 
mechanism found in [27, 28]. In this process, an MHD 
outflow induced by the magneto-spring mechanism is first 
driven along the rotation axis. This outflow then triggers 
convective motions, leading to the formation of large- 
scale eddies in the meridional plane which wind up the 
poloidal magnetic field lines. The strengthened poloidal 
field then leads to further amplification of the toroidal 
field. In our simulations, however, we do not find notice- 
able convection. The plausible reason for this discrep- 
ancy is that the simulations of [27, 28] use a treatment of 
the equation of state and microphysics, which differs from 
ours. In their simulations, neutrino cooling is taken into 
account, and this cooling leads to negative entropy gradi- 
ents and subsequent convection near the neutrinosphere. 
In order to see this convective- type MRI [16], we would 
have to take neutrino processes into account. This is an 
issue to be investigated in the future. 

The simulations presented here were carried out as- 
suming equatorial plane symmetry, and the center of the 
PNS thus remains at the origin. In the absence of this 
symmetry, the PNS may move, due to back-reaction of 
MHD outflows. The lack of equatorial symmetry is cru- 
cial in the explosion mechanism proposed by [2] , in which 
acoustic waves from I = 1 <?-modes provide an impor- 
tant source of energy. In addition, the outflow may de- 
velop anisotropically, since it is driven in part by inhomo- 
geneous MRI-induced turbulence. Anisotropic outflow 
could also arise if the magnetic field profile is anisotropic 
in the supernova progenitor. Our numerical results show 
that <; O.IMq of material will be ejected in the first 
~ 100 ms after the growth of the toroidal field saturates. 
The typical outflow velocity is <~ 0.1c. If the anisotropy 
in the direction of the ejected mass is <; 10%, the PNS 
will move with a velocity <; 10~ 3 c ~ 300 km/s due to 
the back-reaction. Anisotropic MHD outflows may thus 
be able to explain high- velocity pulsars [82] . Simulations 
without equatorial plane symmetry would be required to 
explore this possibility. 

An additional limitation of our simulations is the as- 
sumption of axisymmetry. Nonaxisymmetric instabili- 
ties (such as bar modes and/or one-armed spirals) may 
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arise in the formation of PNSs and contribute strongly 
to the gravitational wave signal [83, 84]. In this paper, 
we mainly consider the case in which the progenitor is 
rigidly rotating and hence the resulting PNSs are weakly 
differentially rotating with T Iot /\W\ < 0.12. This im- 
plies that the PNSs found in this paper would be sta- 
ble against nonaxisymmetric deformations. However, for 
more rapidly rotating PNSs, the possibility of nonax- 
isymmetric deformations have to be taken into account. 
The behavior of the MRI is also expected to be different 
in a full 3D calculation because of the effect of nonax- 
isymmetric MRI induced by toroidal magnetic field [85]. 
Turbulence may arise and persist more readily in 3D 
due to the lack of symmetry. More specifically, accord- 
ing to the axisymmetric anti-dynamo theorem [86], sus- 
tained growth of the magnetic field energy is not possible 
through axisymmetric turbulence. Simulations in full 3D 
will eventually be necessary in order to fully understand 
the role of the MRI in PNS evolution and jet formation. 
Given current computational resources, however, we con- 
sider this a challenge for the future. 
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APPENDIX A: MULTIPLE TRANSITION 
FISHEYE COORDINATES 

The multiple transition fisheye coordinates [34] x l are 
related to the original coordinates x % through the follow- 
ing transformation: 



x <-\ 



r(r) 



a„ r - 



5^K» In 



cosh[(r + r 0i )/sj 
cosh[(f - f 0i )/si 



(oi-i - ai)Si 
2tanh(f 0l /s l ) ' 



(Al) 



(A2) 



(A3) 



where r = y/x 2 +y 2 + z 2 , f — \Jx 2 + y 2 + z 2 , n, ai, roi 
and Si are constant parameters. We perform simulations 
on this fisheye grid with < x < Z, and < z < L. The 
Cartoon method is used to impose axisymmetry as usual. 
We use a grid size of N x 3 x N with uniform grid spacing 
A = L/N. It follows from the transformation (Al) that 
this corresponds to a resolution in the original grid of 
A w (dr/df)A rts a^A in regions well separated from 
the transitions, i.e. where foi + s, <C f <C foi+i — s»+i. 



The ratio A/ A smoothly changes from a% to dj+i in the 
transition region f 0l+1 - s i+ i < r < f i+i + s i+ i. 

Our implementation is as follows. We first perform 
simulations in the original coordinates x l with the re- 
gridding technique as discussed in Sec. V until the value 
of <f> c reaches 0.16. At this time, we interpolate the data 
to the fisheye coordinates x l with the fisheye parameters 
n = 3, (a ,ai,a 2 ,a 3 ) = (0.1,0.4,0.8,1.0). We choose 
the value of L so that the outer boundary is approxi- 
mately the same as that of the original grid. We set 
N = 600, and the values of f oi so that (r i, r 02 , r 03 ) w 
(100 km, 200 km,400 km), where r j = r(f oi ). The 
transition width is set to be Sj = 29A. With this set- 
ting, the resulting resolution in the original grid becomes 
A w 0.7 km for r < 100 km, A w 2.7 km for 100 km 
< r < 200 km, A w 5.5 km for 200 km < r < 400 km, 
and A w 6.8 km for r > 400 km. Therefore, using this 
technique, we can achieve high resolution in the central 
region with relatively few grid points. 



APPENDIX B: PROOF OF k = p,v p /B p 
CONSTANT ALONG MAGNETIC FIELD LINES 

In a stationary spacetime, d t B l = = d t p*. It fol- 
lows from the induction equation (31) and the continuity 
equation (32) that 



d j (v>&-v i &) = 0, 
dj{p*v j ) = 0. 



(Bl) 
(B2) 



In an axisymmetric spacetime, the continuity equation 
becomes 



dp{mp*v p ) = 0, 



(B3) 



where P denotes the poloidal components {w and z). We 
introduce a one-form 



w» = e llk v j B k = [ijk]v j B k , 



(B4) 



where e a( g 7 = n M e MQ( 3 7 is the 3-dimensional Levi-Civita 
tensor, and [ijk] is the permutation symbol. It is easy 
to show that Eq. (Bl) is equivalent to D^uu^ — (i.e. 
uji is a closed one- form), where Di denotes the covariant 
derivative associated with the three-metric ■jij . Hence u>i 
can be written as 



Dif = dif, 



(B5) 



where / is a scalar function. We assume that Eq. (B5) 
holds globally. It follows that uj v = m(v z B™-v™B z ) = 
in an axisymmetric and stationary spacetime. This im- 
plies that /B™ — v z /B z , i.e. v p = fiB p , where fi is 
a scalar function. Substituting v p = \iB p into Eq. (B3) 
we have dp(mp*fiB p ) = 0. Using the no-monopole con- 
straint dp(vjB p ) = 0, we obtain B p dpk = = B^djk, 
where k = = p*v p /B p . Hence k is constant along 
magnetic field lines in an axisymmetric, stationary space- 
time. 
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